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Abstract 

Cleveland State University (CSU), the University of Minnesota (UMN), Gedeon Associates, Infinia 
Corp., Sunpower Inc. and International Mezzo Technologies (Mezzo) have completed Phase II of a 
Radioisotope-Power-Conversion Technology, NASA Research Award (NRA) contract. The project 
brought together experts in Stirling-cycle machine design, microfabrication processing, oscillatory-fluid- 
flow experimentation and Computational Fluid Dynamics (CFD) to design and fabricate an advanced 
regenerator matrix for use in Stirling-cycle space-power conversion technology. The main objectives 
were to significantly increase the overall thermal efficiency of the regenerator and the Stirling convertor 
and to improve the structural reliability and manufacturability of the regenerator. 

Regenerators for the next generation of Stirling convertors should have microscale features that have 
been configured for improved reliability, heat transfer and fluid flow characteristics. These microscale 
features can be produced by batch-mode Electrical-Discharge Machining (EDM) and LiGA (X-Ray 
Lithography) processes. These processes utilize technologies that have been developed for 
MicroElectroMechanical Systems (MEMS). 

During Phase II an actual-size microfabricated regenerator comprised of a stack of 42 disks, 19 mm 
diameter and 0.25 mm thick, with layers of microscopic, segmented, involute-shaped flow channels was 
fabricated and tested. The geometry resembles layers of uniformly-spaced segmented-parallel-plates, 
except the plates are curved. Each disk was made from electro-plated nickel using the LiGA process. For 
historical reasons the geometry is sometimes referred to as “involute-foil.” This regenerator had feature 
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sizes close to those required for an actual Stirling engine but the overall regenerator dimensions were 
sized for the NASA/Sunpower oscillating- flow regenerator test rig. Examination by scanning electron 
microscope showed that the disks were an accurate rendition of the design specification except for a few 
flaws which were expected to be worked out of the manufacturing process. Testing in the oscillating-flow 
test rig showed the regenerator performed extremely well, producing the highest figures of merit ever 
recorded for any regenerator tested in that rig over its ~20 years of use. 

Progress was also made in understanding the detailed fluid dynamics in the regenerator by CFD 
analysis at Cleveland State University and large-scale testing at the University of Minnesota. In general, 
the conclusions of CFD and large-scale testing reinforced the actual-size test results and revealed some 
important details about the microscopic flows and heat transfer responsible for the overall regenerator 
behavior. 

The original batch-EDM regenerator fabrication process anticipated in Phase I did not appear to be 
fast or controlled enough to yield a high quality product in a reasonable time. Therefore, fabrication by a 
FiGA electro-plating process was chosen for the Phase II nickel prototype test regenerator; however 
higher-temperature lower-conductivity more-robust materials, than pure nickel, are required. 
Electroplating from low conductivity and corrosion-resistant materials (even alloys) is possible, but has 
not been explored under this effort. 

In Phase III (now underway) the plan is to microfabricate and test an involute-foil regenerator in a 
Stirling engine. The two key elements of this work are to adapt a Sunpower space-power engine to use an 
involute-foil regenerator, and to fabricate and test the regenerator in the engine. This regenerator will also 
be an electroplated nickel structure. 

For possible work beyond Phase III, see the conclusions and recommendations at the end of this 
report. 


1.0 Introduction 

The Stirling-engine regenerator has been called “the crucial component,” Organ (2000), in the 
Stirling-cycle engine. The regenerator, which obtains heat from the hot working fluid and releases heat to 
the cold working fluid, recycles energy internally, allowing the Stirling cycle to achieve high efficiency. 
The location of the regenerator within a Stirling convertor is shown in figure 1.1. 

Currently, regenerators are usually made of woven screens or random fibers. Woven-screen 
regenerators have relatively high flow friction. They also require long assembly times which tends to 
increase their cost. Random fiber regenerators also have high flow friction but are easy to fabricate and 
therefore are inexpensive. Figure 1.2 shows a typical random-fiber regenerator and figure 1.3 shows a 
close up of the fibers. Due to the method of fabrication, the fibers are random primarily in a plane 
perpendicular to the main flow path. Thus both woven screens and random fibers experience flow 
primarily across the wires (cylinders in cross flow). Cylinders in cross flow tend to cause flow separation 
resulting in high flow friction and considerable thermal dispersion, a thermal loss mechanism that causes 
an increase in apparent axial thermal conduction. For space engines, there must be assurance that no 
fibers of this matrix will eventually work loose and damage vital convertor parts during the mission. It is 
also important that local variations in porosity inherent to random fiber regenerators not result in local 
mismatches in flow channels which would contribute to axial thermal transport. Wire screens have some 
randomness associated with their stacking and thus may have locally non-uniform flow. The efforts thus 
far have shown that attractive features for effecting high fluid-to-matrix heat transfer with low pressure 
drop are a matrix in which: a) the heat transfer surface is smooth, b) the flow acceleration rates are 
controlled, c) flow separation is minimized and d) passages are provided to allow radial mass flow for a 
more uniform distribution when the inlet flow or the in-channel characteristics are not radially uniform. It 
is thought that properly designed microfabricated regular geometries could not only reduce pressure drop, 
maintain high heat transfer and allow some flow redistribution when needed, but could show improved 
regenerator durability for long missions. 
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Figure 1.1 . — Schematic of Stirling convertor showing the location of the regenerator. 




Figure 1.2. — Random Fiber Regenerator. 



Figure 1.3. — Electron micrography of a random fiber regenerator matrix. Courtesy of NASA Glenn Research Center. 
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In Phase I of this project, a microscale regenerator design was developed based on state-of-the-art 
analytic and computational tools. For this design, a 6 to 9% engine-efficiency improvement was 
projected. A manufacturing process was identified and a vendor (International Mezzo Technologies) was 
selected to apply it. Mezzo completed EDM tools for fabricating layers of the chosen involute-foil 
microregenerator design, based on the team’s specifications. They were ready to begin producing 
regenerator layers (annular portions of disks). Also, a Large-Scale-Mock-Up (LSMU) involute-foil 
regenerator was designed and fabrication had begun. Computational Fluid Dynamics (CFD) for different 
geometries was employed to model the fluid flow and heat transfer under both steady and oscillatory-flow 
conditions. The effects of surface roughness were included. Several geometries: lenticular, parallel plates 
(equally/non-equally spaced), staggered parallel plates (equally /non-equally spaced) and 3-D involute- 
foils were studied via CFD. The modeling was applied to both the microscale involute-foil regenerator 
and to the LSMU model of it. 

A number of action items were addressed at the first-year review meeting held at NASA Glenn 
Research Center (GRC) on August 27, 2004. These items and their proposed resolutions were addressed 
in the final report for Phase I. One of these action items requested clarification of a plan for structural 
analysis of the microfabricated regenerator. Preliminary analysis showed that the involute-foil regenerator 
is of a sound design, structurally. Extremely high axial stiffness ensures that with an appropriate axial- 
compression force, the regenerator can be held firmly in place with negligible regenerator deformation. 
Special caution is needed for regenerator installation to prevent potential lateral deformation due to 
misalignment, as the radial stiffness is relatively low. 

The goal of the current NASA project is to develop a new regenerator of high durability as well as 
high efficiency using emerging microfabrication technology. In addition to the benefit to Stirling 
convertor space-power technology, such regenerator development will also benefit Stirling cycle coolers 
and NASA’s many cryocooler-enabled missions. 

This report will show how the microscale involute-foil regenerator was fabricated and tested, and will 
review the test results. Also, the LSMU experimental model and the data obtained from it will be discussed. 
Moreover, the newly designed regenerator was analyzed via 1-D, 2-D and 3-D computational analyses . 

These computational tools enabled evaluating: 1) The figure of merit and 2) The effect of varying different 
parameters such as solid material used, thermal contact resistance, etc. for the new design. 

2.0 Nomenclature 

A, Aj Area, m 2 

b Jet width, m 

C Specific heat, J/kg.K 

CA Crank Angle, degrees 

Cf Inertial coefficient 

D/„ du Matrix hydraulic diameter, m 

D p Piston bore diameter, m 

/ Oscillating frequency, Hz, or Darcy friction factor 

F Fraction of matrix not participating fully in thermal exchange with fluid, but equation 

(3.9) defines another usage of “F” 
h Convective heat transfer coefficient, W/'nT K 

K Permeability, m 2 

K, kj' Thermal conductivity, W/m.K 

L Length, m 

Nu Nusselt number 

Nu_m Mean Nusselt number 

Nu x Local Nusslet number 

N k “conductivity ratio” defined as the effective axial conductivity divided by the molecular 

conductivity 
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Re * Pr = Peclet number 
Prandtl number 
Pleat into engine, W 
heat transfer rate, W/m 2 
Pressure drop. Pa 

Radial distance from the test section axis, m 

Radius, m 

Reynolds number 

Maximum Reynolds number 

Spacing between two cooler tubes, m 

time, s 

Temperature, °C 
Cold end temperature, °C 
Plot end temperature, °C 

Maximum bulk mean velocity in regenerator, m/s 

Darcy velocity in streamwise direction, m/s 

Volume, m 3 , or local velocity inside the channels, m/s 

PV, or indicated, power predicted for engine, W 

Displacer drive power, W 

Valensi number 

Axial distance, m 

The jet penetration depth, m 

Amplitude of piston displacement, m 

Amplitude of particle displacement within jet generator tubes, m 
Amplitude of particle displacement within regenerator, m 
Maximum particle displacement amplitude, m 
Porosity 

Crank angle, degrees 
Dynamic viscosity, kg/m sec 
Kinematic viscosity, m 2 /sec 
Density, kg/m 3 

Stefan-Boltzman Constant, W/m 2 -K 4 
Dimensionless temperature or porosity of porous media 
Angular frequency, rad/sec 

2,1 Abbreviations 

Advanced Stirling Convertor (engine and linear alternator) 

Computer Aided Design 
Computational Fluid Dynamics 
Cleveland State University 
Electric Discharge Machining 
Finite Element Analysis 

Frequency Test Bed (convertor, or engine and linear alternator) 

Glenn Research Center 

Lithographie, Galvanoformung and Abformung (the German words for lithography, 
electroplating and molding. X-ray lithography is used here) 

Large Scale Mock Up (of involute-foils) 
microelectromechanical systems 
National Aeronautics and Space Administration 
NASA Research Award 
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PMMA polymethyl methacrylate (a clear plastic, also marketed as Acrylic, Plexiglas, Lucite, etc. 

Used as a photoresist in LiGA process for microfabrication of involute-foils) 

SEM Scanning Electron Microscope 

SU-8 SU-8 is a negative, epoxy-type, near-ultraviolet photoresist used in 

microelectromechanical (MEMS) applications 
TCR Thermal Contact Resistance 

UMN University of Minnesota 


3.0 Program Overview 

3.1 Objectives 

The main objective of this phase (II) was to build a prototype of the involute-foil regenerator and test 
it in a NASA/Sunpower oscillatory flow rig to determine how well it performs relative to currently used 
regenerator structures like random fiber and wire screen. This was accomplished using modem 
microfabrication techniques supported by large-scale involute-foil testing, theoretical analysis and 
computational-fluid-dynamics (CFD) modeling. 

Forty-two layers of involute-foil disks (i.e., annular portions of disks) were microfabricated via a 
LiGA process; LiGA is derived from the german words for lithography, electroplating and molding. The 
large-scale testing was based on dimensions 30-times actual size, and test conditions dynamically similar 
to those expected in an engine. The CFD computations were based upon 2-D and 3-D simulations of 
representative elementary regenerator volumes, and 1-D simulations of regenerators and engines. For the 
Stirling-space-engine application, goals for the new regenerator were improved performance, reliability 
and manufacturability. 


3.2 Approach 

The approach included: 1) reviewing the literature, 2) identifying and computing performance figures 
of merit, 3) utilizing state-of-the-art CFD modeling and large-scale involute-foil testing, 4) applying state- 
of-the-art microfabrication techniques and 5) experimentally testing the fabricated microscale regenerator. 
A brief description of these items is given below. 

First, a review of published theoretical and experimental regenerator performance literature was 
completed. Computational research was done in the area of microchannel fluid-dynamics and heat 
transfer to learn about the effect of roughness in small channels. In addition, a review of different 
manufacturing techniques for micro channels was carried out. Performance figures of merit were 
identified, and improved upon, based upon analytical solutions for simple channel geometries as well as 
from experimental test data available in the literature. State-of-the-art modeling was utilized via a 1-D 
Sage Stirling cycle simulation model and 2-D and 3-D CFD Navier-Stokes solvers for simulation of 
microscale regenerator geometry. In the microfabrication area, available methods were surveyed to 
provide the required matrix within the planned time and cost. Design, fabrication and experimental testing 
of a LSMU (Large Scale Mock Up) which was dynamically similar to an engine involute-foil matrix was 
carried out. Also, testing of a microscale involute-foil was done in an existing NASA/Sunpower- 
developed oscillating-flow test rig. 


3.3 Technical Challenges 

The technical challenges faced in this project are summarized below: 

1) Microfabrication of Small-Feature Sizes: These features are sub- 100 pm scale that must be built in 
an overall macro-scale regenerator of 100 mm scale. No more that a 10% passage-width nonuniformity 
tolerance (10 pm) in these channels was specified, based on computations made by Gedeon Associates. In 
addition, a high passage aspect ratio (length > 3 x width) for each element of the regenerator was 
specified, which required some research to achieve. If the originally planned EDM had been used, in 
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addition to LiGA, EDM tool wear and positioning would be limiting factors. Cost was also limiting, due 
to the limited funding. 

2) CFD of Complex Geometries: There were computational challenges in dealing with the actual 3-D 
overall structures having interrupted passages. Examining the effects of surface roughness (on heat 
transfer) in micro channels required fine resolution of near- wall flow and heat transfer. These challenges 
are demonstrated by the resulting requirements of (a) the large number of computational grid points (for 
proper accuracy of the solution) and, (b) the long computation times, particularly for the oscillatory-flow 
simulations. 


3.4 Phase II Work Plan Review 

The following tasks were carried out during Phase II: 

Task 1 — Large Scale Mock Up (LSMU) final design, test plan and fabrication. 

Task 2 — Testing of LSMU. 

Task 3 — Continued computational modeling of LSMU and microscale involute-foils. 

Task 4 — Fabrication of electroplated nickel prototype regenerator matrix with actual-size features. 

Task 5 — Scanning Electron Microscope (SEM) inspection of microscale involute-foils. 

Task 6 — Actual-size oscillatory-flow testing in the NASA/Sunpower Test Rig. 

Task 7 — Documentation via technical and financial reporting, oral review presentations and a final 
report. 

The following milestones were established to track Phase II progress: 

MS 1 : Finalize design of the LSMU. 

MS2: Complete fabrication of the LSMU and complete the LSMU test plan. 

MS3: Complete fabrication of actual-size prototype regenerator and assessment of resulting matrix. 

MS4: Complete initial testing of the LSMU. 

MS5: Testing of the microfabricated nickel regenerator. 

3.5 Description of Technology 

3.5.1 Regenerator Concepts 

In Phase I, different concepts were identified for the regenerator matrix that could be manufactured 
with existing technology and be expected to have better performance than existing matrices. These 
concepts included: 1) a “lenticular” array, 2) a honeycomb structure, 3) an involute-foil and 4) a modified 
involute-foil. 

The modified involute-foil was selected. This design has multiple channels of uniform width with 
good flow and heat transfer characteristics. It is durable. In contrast to the foil (or wrapped-foil) 
regenerator, although the chosen regenerator concept does have some of the low-loss attributes of the foil 
regenerator, it has uniformly repeating geometric features with a robust structure that prevents the 
variability in channel size that has proven inherent in the wrapped-foil regenerator. In comparison to the 
random-wire regenerator or the screen regenerator, the chosen regenerator concept has the durability of 
the screen regenerator but with better flow-loss characteristics. This modified involute-foil regenerator 
can be fabricated with modem microfabrication techniques. 

3.5.2 Selected Regenerator Concept 

This regenerator consists of a stack of involute-foil disks (or annular portions of disks) that have been 
microfabricated by the LiGA or LiGA/EDM processes (the particular material determines whether EDM is 
required). It has a low resistance to flow because it has a reduced number of separation regions, compared to 
wire screen and random fiber. But it still has high heat transfer effectiveness approaching that of the other 
structures. The resulting figure of merit (~ the ratio of heat transfer to pressure drop) has proven superior to 
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the currently used random- fiber and wire-screen regenerator structures. It has known small features which are 
secured to the assembly so that it is not as susceptible to releasing small fragments of regenerator matrix 
material (as random- fiber and wire-screen, both of which have random features). 

3.6 Summary of Phase II Accomplishments 

3.6.1 Prototype Regenerator Fabrication (Mezzo) 

The LiGA micromachining process was used to fabricate a nickel involute-foil regenerator that was 
tested and found to provide very good performance. The original manufacturing approach of using LiGA- 
fabricated EDM tools to fabricate stainless-steel regenerator parts seemed initially to offer little potential 
due to the extremely low material removal rate. In retrospect however, if higher-machine bum-rate 
settings and EDM tool dimensions are chosen such that greater overbum is allowable, then the EDM- 
LiGA approach will be more viable. It should also be noted that, in this effort, a new manufacturing 
technique was developed: namely using EDM to “skim cut” regenerator parts fabricated by LiGA. While 
the regenerator tested did provide excellent performance, LiGA and/or a LiGA-EDM process 
optimization could result in a better product. Potential improvements include: 

i) Improve the lithography process to eliminate or greatly reduce “undercutting.” 

ii) Cease electroplating before overplating begins. This would eliminate the need to use the “skim 
cut” and would eliminate the source of burrs that are attributed to the EDM operation. 

iii) To greatly reduce cost, explore the use of SU-8, a negative resist that requires substantially less 
exposure time than PMMA. 

iv) For EDM, find acceptable combination of material removal rate, overbum, and geometry that 
gives a high quality part in a reasonable time and a reasonable cost. 

3.6.2 Prototype Regenerator Testing (Gedeon Associates/Sunpower) 

3.6.2.1 Physical Description 

The regenerator sample tested consisted of a stack of 42 involute-foil disks (layers). The following CAD 
rendering (fig. 3.1) shows progressive magnifications of a typical disk viewed from the front. 

The term “involute” refers to the curved shape of the channels. The centerline of each channel might be 
imagined as traced by a fixed point on a string unwinding from a circular cylinder (the generating circle) 
concentric with the disk. The diameter of that cylinder is equal to the diameter of the inner circular rib for the 
previous ring of flow channels. 

The next CAD rendering (fig. 3.2) shows a typical single flow channel in the matrix. It is from an early 
solid model and does not correspond exactly to the final matrix geometry but it is useful to illustrate some of 
the important dimensions that define the geometry. L, is the flow channel length (disk thickness), W is the 
channel width (arc length of web), g is the channel flow gap (normal distance across) and 5 is the basic 
involute element spacing (gap + web thickness). 

From the CAD drawings it is possible to calculate the representative hydraulic diameter for the entire 
regenerator matrix like this: A unit cell of the matrix consists of one 7-ring disk and one 8-ring disk. The 
hydraulic diameter for such a 2-disk cell is a reasonable approximation for the entire matrix. There are 7 
different basic channel shapes in the 7-ring disk and 8 in the 8-ring disk. For each basic channel shape we use 
available CAD tools to measure individual flow area and wetted perimeter. The total flow area A T for the 2- 
disk cell is the sum of the areas of each individual channel elements multiplied by the number of occurrences 
per disk. The same is done to compute the total wetted perimeter W T . The number of channel element 
occurrences per disk is just its involute generating circle diameter divided by the 100 pm involute spacing^ 
(see app. D). The final representative hydraulic diameter is 4A t /Wt .The final hydraulic diameter 
calculated this way is: 


D h = 162.0 pm 
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Figure 3.1. — Involute-foil geometry. Figure 3.2. — Channel geometry. 


The above value agrees within 0.2% of the value computed by UMN for the large-scale mock-up 
(4.872 mm) when scaled up by the scale factor 30. For parallel plates, the hydraulic diameter would be 
exactly twice the flow gap ( 2g) which is more like 170 pm. 

For overall matrix porosity, the value reported by Mezzo was used: 

p = 0.8384 

This value was based on physical measurements of mass and the known material density for nickel. 
As reported by Mezzo, it agrees quite well with the theoretical average porosity for a 2-disk cell of 
0.8299. There is only a 1% porosity discrepancy in the direction that suggests that the flow channel gaps 
are slightly wider than expected (by about 1%). From this observation, one might expect the actual 
hydraulic diameter to be larger by about the same amount, on average, which would bring it to about 
163.7 pm (Di, scales with flow area if wetted perimeter remains about the same). 

Another parameter of interest is the ratio of channel flow length L c (disk thickness) to hydraulic 
diameter which CSU has been investigating with CFD modeling (see sec. 6.0). For the current batch of 
disks, the mean thickness is 238 pm (42 disks totaling 10.0 mm ) so: 

L c ID h = 1.47 

Still another parameter of some importance is the channel aspect ratio, or ratio of channel width to 
hydraulic diameter. The weighted average flow channel width for the two types of disks is about 1200 pm 
so the average channel aspect ratio is: 


W/ D h = 7.4 

In addition to the above parameters, the way the disks are stacked on top of each other is important. 
The current scheme of alternating involute directional orientation and staggered ring walls is easy to 
describe but difficult to quantify in terms of any simple numerical ratio. 

3. 6.2.2 Test Canister Design 

The regenerator disks are housed within the test canister depicted in figure 3.3 for testing in the 
NASA/Sunpower oscillating-flow rig: 
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involute disk stak 



0.8 UMN flat head screws) 

Figure 3.3. — Testing canister geometry. 

3.6.2.3 Regenerator Verification from SEM Images 

Scanning electron microscope (SEM) images of the involute-foil regenerator disks revealed excellent 
spacing uniformity of the flow channels and a generally smooth surface finish but also some defects. 
SEM micrographs were taken with the assistance of Randy Bowman at NASA GRC. 

3. 6.2.3. 1 The Specimen 

The sample photographed was the first one tested in the NASA/Sunpower oscillating-flow test rig. 
The matrix was not taken apart. The purpose was not to do an exhaustive survey of every one of the 42 
disks inside but only to take a look at what was visible from the two ends of the assembled regenerator 
canister. 


3.6.2.3.2 Spacing Uniformity 

One objective was to assess the uniformity of spacing of the foil elements in the matrix. It was 
previously established that the spacing uniformity should be within ±10% to avoid significant adverse 
impact on the figure of merit (app. B). It was found that the Mezzo matrix is significantly better than this, 
as can be seen in the following images. The first image (fig. 3.4a) is a low-magnification view that gives 
one the qualitative impression that the spacing is uniform. The second image (fig. 3.4b) shows actual 
measurements with the Adobe Acrobat “measurement” tool used to conclude that the spacing uniformity 
is within ±2% in a small region of the matrix. The spacing may actually be more uniform than this 
because there was a significant error in the accuracy with which the estimate of the normal direction 
across the channel and the location of the channel edges was done with the measurement crosshairs. 

The measurements are hard to read on the photograph but the values were saved in a “csv” (comma- 
delimited text) file (option in measurement dialog). Using Microsoft Excel’s (Microsoft Coiporation) 
statistical functions the standard deviation of foil spacing was computed to be 1.6% of the average 
spacing, for the 12 measurements. 
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Figure 3.4. — (a) Face-on view at low magnification showing overall good foil 
spacing uniformity. 



Figure 3.4. — (b) Adobe Acrobat measurement tool in action measuring localized foil 
spacing in a close-up image. 
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3.6.2.3.3 Perspective View 

The images below were obtained by tilting the sample by 25°. The first one (fig. 3.5) gives a good 
idea of the 3-D matrix structure and also shows a very smooth surface finish on the walls of the flow 
channels. The second (fig. 3.6) shows a magnified view of a channel wall. We estimated the “roughness” 
height to be less than 1 pm, compared to an 85 pm channel gap. The roughness mainly consists of 
occasional pores and shallow grooves or steps, parallel to the flow direction. 



Figure 3.5. — A 25° tilted view showing 3-D structure near an annular ring with 
involute webs above and below. The dust visible on the surface is 
probably from testing or handling. 



Figure 3.6. — Detail of flow channel surface roughness. 
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Figure 3.7. — Visible debris attached to the face of the second disk within the sample, 
visible below the top disk. 

3.6.2.3.4 Defects 

The matrix was not without its flaws. Figure 3.7 shows significant spatter-like debris on the face of the 
second disk below the surface of the regenerator — at one location (Microscope operators tend to focus on the 
defects). The debris comes from EDM removal of the disk from the plating substrate. The top disk also shows 
some evidence of this debris on its lower surface. Other images show similar debris comprised of small 
spherical particles of what appears to be nickel (preliminary SEM assay). 

Eliminating the debris and edge roughness would be an important objective in further development. 
Roughness will have some effect (documented in the oscillatory flow tests) on the helium entering the 
flow channels and might pose a contamination problem if any of the material works loose. 

Another “flaw” in this particular assembly is that the angular involute orientations are the same in the 
top two disks instead of opposite (crossed), as they are supposed to be. Subsequent inspection by Mezzo 
showed that there was a more-or-less random involute orientation throughout the matrix, although the two 
types of disks were correctly alternated. This stacking problem was corrected and the properly-stacked 
regenerator was also tested. 

There is another sort of defect (see fig. 3.8) that appears to be associated with a flaw in the photo- 
lithographic mask used to expose the photo-resist material. Occasionally, one sees notches extending the 
full disk thickness. This particular notch extends a bit less than halfway through the web thickness. The 
presence of such defects might put an effective lower limit on the thinness of the webs before a significant 
number are cut completely through. In fact, it may explain why some webs were cut completely through 
in some cases. 

3.6.2.3.5 Recommendations 

In spite of the flaws, the test results for this matrix were very encouraging. However, some of the 
flaws are of concern for long-term reliability so work on quality control issues is needed for future 
matrices. Based on these photographs the two main concerns are spatter-like debris near the ends of the 
flow channels and notch defects extending completely through or mostly through channel walls. 
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Figure 3.8. — Notch defect extending for the full length of an involute web. 


Figures of merit — various matrices 



Reynolds number 

Figure 3.9. — Figures of merit for various matrices. 


3.6.2.4 Test Results 

The results of testing in the NASA/Sunpower oscillating flow test rig were very promising and are 
summarized along with results from several other regenerator types in figure 3.9. 

The microfabricated regenerator has a figure of merit substantially higher than the other regenerator 
types — including the 90% random fiber regenerator which is roughly what is being used in the current 
generation of space-power Stirling engines. Figure of merit is defined as follows: 


NASA/CR— 2007-2 15006 


14 


1 


R'Pr , N t ] 

4^„ W 

■R, Reynolds number 
Prandtl number 
f Friction factor 

N u Nusselt number 

N k Thermal dispersion conduction 
enhancement 


M / 

f 

V 


The figure of merit is a first-cut measure of overall regenerator performance. It is inversely 
proportional to the product of regenerator pumping loss, W p , thermal loss, Q t , and the square of 
regenerator mean flow area, A f (see details in app. F, by Gedeon), as follows:. 


F m <*■ 


W p Q t Aj 


T/ tends to be constrained by power density (void volume) so when comparing regenerators in similar 
engines it may be ignored. 


3.6.2.4.1 New 96% Porosity Random Fiber Data 

In the plot (fig. 3.9) the figure of merit for 96% porosity random fibers is based on recent test data 
with improved accuracy compared to data previously reported. Based on the earlier data, it appeared that 
the figure of merit for 96% porosity random fibers was actually higher than that of the microfabricated 
regenerator at some Reynolds numbers. That is no longer the case and the microfabricated regenerator 
now has a figure of merit substantially higher than that for the 96% porosity random fibers. The 
microfabricated regenerator now ranks as the best regenerator ever tested in the NASA/Sunpower test rig. 


3.6.2.4.2 Unexplored Dimensionless Ratios 

The pressure-drop and heat transfer correlations discussed below (eqs. (3.1) through (3.4)) are for a 
particular regenerator matrix. Since the correlations are in terms of dimensionless quantities (like 
Reynolds number) they also apply to geometrically similar matrices. What does that mean? 

Aside from the essential geometric property of the involute channels, namely that they consist of 
uniform-gap planar flow passages, and the chosen stacking geometry, the important dimensionless 
specifications are porosity p, the ratio of flow channel length to hydraulic diameter, LJD h , and the aspect 
ratio, W/D ) Of these probably LJD h is most important because it affects the degree to which flow is fully 
developed within any given flow channel at a given Reynolds number. Porosity is probably of lesser 
importance because it does not directly affect the nature of the flow channels. But it does affect the 
thickness of the webs between flow channels and therefore the way the flow enters the channels in 
distributing from one disk to the next. The aspect ratio is probably not too critical so long it is not 
significantly lower than the value mentioned above. 

So one should beware of applying the correlations below to microfabricated involutes with 
significantly different porosity, LJD h or W/D h . But if one does, it is expected that the present correlations 
are conservative when applied to regenerators with higher porosity or higher L c !D h or higher W/D h In 


'Since the porosity is a relatively good measure of the ratio of flow gap to element spacing, g/s, and the fill factor 
( 1— P) is a relatively good measure of the ratio of web thickness to element spacing (l-g/s) there is no need to 
separately use g/s or (1-g/s) to characterize the matrix. 
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either of those cases the flow would then be closer to the ideal parallel plate flow and the figure of merit 
should increase somewhat. 

3.6.2.4.3 Restacked Canister 

Since the original test regenerator was found to have been incorrectly assembled (random spiral 
orientation) after it had been tested in the oscillating-flow test rig, it was decided that the spiral orientation 
was important enough to warrant re-testing the regenerator. Mezzo re-stacked the regenerator using the 
same disks but this time with disks both correctly sequenced and with the spiral direction reversed at each 
disk transition, according to the original plan. It was then tested in the oscillating-flow rig and found that 
the overall figure of merit changed slightly compared to the original testing, as shown in figure 3.10. 

Figure 3.10 shows that the correct stacking produces a slightly better figure of merit at high Reynolds 
numbers and slightly worse values at low Reynolds numbers. As will be seen below, the friction factors 
for the two cases are almost identical so the main reasons for the differences are thermal losses, with the 
restacked regenerator producing more thermal loss at low Reynolds numbers and less at high Reynolds 
numbers. The increased thermal loss at low Reynolds numbers is probably a result of improved 
measurement accuracy as a result of new rig operating procedures, as explained below. The reduced 
thermal loss at high Reynolds numbers appears to be real. 

To get a better feel for what is going on two data points near the peak of the figure of merit curve, at 
Reynolds numbers around 400 were taken, one data point for the original regenerator and one for the re- 
stacked canister. Both correspond to tests with 50 bar nitrogen. 

As shown in table 3.1, the two data points are nearly identical in all respects except that the restacked 
canister has about 2.4% lower thermal loss (heat rejection to cooler) for a 2.0% larger regenerator 
temperature difference. Assuming thermal loss scales in proportion to temperature difference, this implies 
that the restacked regenerator would produce about 4.5% lower thermal loss for the same temperature 
difference. This is consistent with the figure of merit calculation which shows about a 5% improvement 
for the restacked regenerator over the originally- stacked regenerator at a Reynolds number of 400. 

Figure of merit — Mezzo Involute Microfab 



Reynolds number 

Figure 3.10. — Figure of merit values, random- and correctly-stacked canister. 


TABLE 3.1— VALUES FOR THE ORIGINALLY- AND CORRECTLY-STACKED REGENERATORS 



Random stacking 

Correct stacking 

Mean pressure (bar) 

50.0 

50.0 

Piston amplitude (mm) 

4.001 

4.000 

Coolant flow rate (g/s) 

6.161 

5.712 

Coolant AT (C) 

2.149 

2.264 

Heat rejection (W) 

55.39 

54.08 

Tmean regen hot end (C) 

449.2 

450.5 

Tmean regen cold end (C) 

340.1 

339.2 

AT regenerator (C) 

109.1 

111.3 
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At low Reynolds numbers, where heat rejection is smaller, the thermal differences between the two 
regenerators are likely due to changes in test-rig operating procedures. For the restacked regenerator a 
“single ramp-up” procedure designed to minimize difficulties with long-term thermal drift was adopted. 
Under this procedure, the operator increases the piston amplitude in 1 mm increments from 0 to 1 0 mm, 
waiting about 30 min for the rig to equilibrate after each change. Under the previous procedure, the 
operator ramped piston amplitude up and down twice over the range 0 to 1 0 mm, with only about 8 min 
equilibration time between changes. There was also a new procedure for measuring the baseline cooler 
heat rejection due to static thermal conduction. The rig was operated at 5 mm piston amplitude (mid 
range) for 2 hr then allowed to sit for 30 min at zero piston amplitude before logging the baseline static 
thermal conduction data points. Previously, static thermal conduction was averaged from zero-amplitude 
data points logged several times during the test without waiting as long for the rig to settle down to 
thennal equilibrium. The new procedure is considered to produce more accurate results at the low 
Reynolds number end of the experimental range. 

3.6.2.4,4 Friction Factor Correlations 

The Darcy friction factors for the original and re-stacked tests are: 

1 20 9 

/ = : — h 0.362 Re -0 056 (original stacking) (3.1) 

Re 

and 

/ _ 117-3 + o 380 Rg— 0 053 ( corr ect stacking) (3.2) 

Re 

Plotted as functions of Reynolds number there is hardly any difference (correct stacking about 0.9% 
lower on average). The following plot (fig. 3.11) focuses in on a small range of Reynolds numbers to 
better resolve the two curves. If it is plotted over the full range from 10 to 1000, then the two curves 
become indistinguishable. 


Darcy friction factor Mezzo Involute Microfab 



Figure 3.1 1 . — Darcy friction factors, f, as functions of Reynolds number, Re, 
for the originally- and correctly-stacked regenerators. 
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The range of key dimensionless groups for these tests was: 


TABLE 3.2.— DIMENSIONLESS GROUPS FOR 


THE PRESSURE DROP TESTS 


Peak Re (Reynolds number) range 

3.4 to 1190 

Va range (Valensi number) 

0.11 to 3.8 

5/L range (tidal amplitude ratio) 

0.13 to 1.3 


3.6.2.4.5 Heat Transfer Correlations: Simultaneous Nu and Nk 

The correlations for simultaneous Nusselt number and enhanced conductivity (thermal dispersion) 
ratio derived for this matrix are: 

Nu = 1 + 1.99Pe 0 358 Nk = 1 + 1.3 14Pe 0358 (random stacking) (3.3) 

and 

Nu = 1 + 1.97Pe 0 374 Nk - 1 + 2.51 9Pe 0 - 374 (correct stacking) (3.4) 

Pe = RePr is the Peclet number. Nu and TV* are plotted individually below for the case Pr = 0.7 (see 
figs. 3.12 and 3.13). 


Mean-Parameter Nu Mezzo Involute Microfab 



Figure 3.12. — Mean-Nusselt-number Nu values as functions of 
Reynolds number, Re, for the originally- and correctly stacked 
regenerators 
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Mean-Parameter Nk Mezzo Involute Microfab 



Figure 3.13. — Mean enhanced-thermal-conductivity (thermal- 
dispersion) ratio, Nk , values as functions of Reynolds number, 
Re, for the originally- and correctly stacked regenerators. 


As a reminder, Nu and N k are designed to be used together in a model like Sage where the Nusselt 
number may be understood as based on section-mean temperature rather than velocity-weighted (bulk) 
temperature. Under this assumption, N k compensates for any discrepancy in enthalpy flow compared to 
using a bulk-temperature approach. 

Comparing results for the two regenerators, in figure 3.12, shows that the Nu values are very close to 
one another except for a slight, but significant, increase at high Reynolds numbers for the correctly- 
stacked case. The N k values in figure 3.13 are quite different, probably due to larger measured thermal 
losses at low Reynolds numbers attributed to lower baseline thermal losses measured under the new rig 
operating procedures. The combined effects of Nu and N k together are best seen in the plot of figure of 
merit at the beginning of section 3. 6.2. 4 (fig. 3.9). 

The range of key dimensionless groups for these tests is: 


TABLE 3.3— TEST PARAMETER RANGES 


Peak Re range (Reynolds number) 

2.6 to 930 

Va range (Valensi number) 

0.064 to 2.4 

5/L range (tidal amplitude ratio) 

0.17 to 1.8 


3.6.2.4.6 Parallel Plate Nusselt Number Comparison 

Figure 3.14 compares two Nusselt-number plots derived for the micro fabricated involutes against the 
theoretical Nusselt number (Nu = 8.23) for fully-developed flow between parallel plates under the 
uniform heat flux boundary condition. 

The curve labeled “Nu” is the Nu part of the simultaneous Nu, N k correlation. The curve labeled 
“Nue” is an effective Nu ( , correlation, derived under the assumption that N k = 1 . The two derived Nusselt 
numbers are close to each other but far from the theoretical Nusselt number. At high Reynolds numbers, 
the derived values are higher than the theoretical value, which makes sense because the flow in the 
microfabricated flow channels is more and more like developing flow with increasing Reynolds numbers 
and Nusselt numbers are known to be higher in developing flow (as documented by CSU see sec. 6.0). At 
low Reynolds numbers the derived values are lower than the theoretical value. This may be due to the 
effects of solid (nickel) thermal conduction within the individual regenerator disks which is of some 
significance at low Reynolds numbers, (see sec. 3.6.5). 
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Nusselt Numbers Mezzo Involute Microfab 



Figure 3.14. — Nusselt-number values as functions of 

Reynolds number, Re, compared to the fully-developed 
channel value 


3.6.3 Large Scale Mockup (University of Minnesota, UMN) 

The UMN contribution in the second year (Phase II) began with the design of a large-scale (30 X 
actual size), dynamically-similar mockup of the microfabricated regenerator for testing with higher spatial 
and temporal resolution than afforded by the actual size regenerator. The geometry, operating conditions 
and instrumentation used in testing with this model are discussed in sections 5.1 and 5.2. The use of 
dynamic similitude was verified by agreement of flow and heat transfer measurements from the Large- 
Scale MockUp (LSMU) with measurements from the NASA/Sunpower oscillatory flow rig. This facility 
was then used to measure frictional pressure drop, time and space-resolved heat transfer rates, and the 
unsteady matrix-flow thermal interaction associated with jets entering the matrix from passages of 
adjacent heat exchangers. Details are given in sections 5.3 through 5.5. 

3.6.3.1 Darcy Friction Factor 

A study of the momentum equation noted that for engine-representative Valensi and Reynolds 
numbers, the transient term is unimportant and pressure drop measurements can be taken in steady, 
unidirectional flow. Such measurements led to the following friction factor correlation for the LSMU: 

/ = — + 0.127 Re 0 01 (3.5) 

Re 

Details are given in section 5.3.3. Figure 3.15 shows the comparison of this correlation to the correlation 
from the NASA/Sunpower oscillatory flow test rig (eq. (3.2)). These two correlations match very well at 
the low end of the Reynolds number range (—100 to 200), with the LSMU correlation no more than ~15% 
higher, and the NASA/Sunpower correlation is about 25% higher at the high end of the Reynolds number 
range (-1000). The LSMU correlation may be higher for lower Re values due to the shortness of the 
LSMU test section. The microfabricated test section friction factor may be higher for larger Re values due 
to the roughness associated with the EDM cutting process, shown earlier in figure 3.7. 

3.6.3.2 Heat Transfer Coefficients 

Unsteady heat transfer measurements were taken in an oscillatory flow with the LSMU experimental 
model. Details of this test are given in section 5.5. The following direct measurements were taken: 1) the 
heat flux between the fluid and the metal matrix was measured for several points interior to the matrix, 
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LSMU plates 

Correlation from 

NASA test rig 


Re 

Figure 3.15. — Comparison of the Darcy friction-factor as functions 
of Reynolds number correlations for the LSMU and the 
microfabricated-actual-size involute-foil regenerators. 


resolved in time within the cycle, and 2) temporally-resolved solid and fluid temperatures were taken at 
the location of the heat flux measurements. Thus, heat-transfer coefficient, h = ” /( n , _ \ and Nusselt 


\Ts~Tf) 


number, Nu = could be computed, resolved in space and location within the oscillation cycle. 

/ K f 


The Nusselt number relationship is complex, as discussed in section 5.5. However, when at and near the 
peak velocity and in the deceleration portion of the cycle, the instantaneous values compare well with 
values computed from a correlation developed from data taken in the NASA/Sunpower oscillatory flow 
experiments with the microfabricated actual-size regenerator (see NASA/Sunpower plot in fig. 5.58, from 
equation (3.4)). The adjusted values in table 5.6 when compared with the microfabricated-regenerator 
Nusselt numbers of figure 5.58, show that the actual-size regenerator values are about 20% higher than 
the adjusted LSMU values. It may be that the microfabricated regenerator values are higher because of 
the roughness shown in figure 3.7, resulting from the EDM skim cutting. Heat transfer rates between the 
gas and the matrix throughout the cycle agree well with values computed from the NASA/Sunpower rig 
correlation, as is shown and discussed in section 5.5 (see fig. 5.59). 


3.6.3.3 Jet Penetration 

One concern with integrating the microfabricated regenerator into the engine is the effectiveness of 
heat transfer on the regenerator ends where discrete jets, formed in the acceptor or rejector heat 
exchanger, penetrate and diffuse within the matrix while exchanging thermal energy with the matrix 
material. To address this, a slot jet and a round jet generator was mated to the LSMU and the thermal 
signatures of those jets were measured within the matrix as the jets (cold in this study) entered into and 
dispersed within the matrix. To supplement the thermal measurements and aid in understanding the 
processes, velocities were measured within the plenum between the LSMU and the jet generator, resolved 
in radial position and in location within the oscillation cycle. Details of this experiment are given in 
section 5.4. No jetting study was conducted in the NASA/Sunpower test rig but CLD computations of the 
LSMU experiments are under way. The computational results will be reported in the Phase III final 
report. 

Two fundamental parameters were extracted from this study, both at the maximum-velocity location 
within the cycle when the jets are immerging into the matrix: 1) the depth into the matrix at which the 
thermal signature of an individual jet can no longer be distinguished (the “jet penetration depth”) and 2) a 
measure of the matrix volume fraction that resides outside of the jets over the matrix volume which 
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extends from the end of the matrix to the penetration depth (the “fraction of inactive material,” F). Since 
the jets diffuse while they penetrate, F is not to be taken too literally. In fact, considerable heat transfer 
between the matrix and the jets occurs beyond the “edge” of the immerging jets. These two values are 
given in table 3.4. 


TABLE 3.4— RESULTS OF THE JET PENETRATION STUDY - LSMU 


Jet geometry 

X / 

Penetration depth, p / , 
/ dh 

(multiples of the matrix 
hydraulic diameter) 

“Inactive” fraction of the matrix 
volume, F 

(between the edge of the matrix 
and the penetration depth) 

Dimensionless total volume of 

X / 

“inactive” matrix, F ' p /, 

/dh 

(nonnalized by the volume Ajd/,) 

Round jet 

13 

0.47 

6.1 

Slot jet 

10 

0.69 

6.9 


A total volume of “inactive” matrix scaled on the volume A jdj, can be computed as: F 



, where 


Aj is the cross-sectional area of the matrix attributable to each jet of the adjoining heat exchanger and 
d /, is the matrix hydraulic diameter. This total volume value is shown also in table 3.4. 


3.6.4 Analysis Tools and CFD Results (CSU) 

The microfabricated involute-foil regenerator was numerically simulated utilizing commercial CFD 
software (Fluent) under both steady-state and oscillatory-flow conditions. The geometry consists of a 
stack of disks with each disk containing involute-shaped micron-range channels, with channel flow 
direction perpendicular to the plane of the disk. The lateral orientation of the channels alternates from 
disk to disk in the flow direction. Simulations were done for both 2-D and 3-D computational domains. 
Steady-state simulations were performed for Reynolds numbers from 50 up to 2000 based on the channel 
hydraulic diameter and the mean flow velocity. Oscillatory flow simulations were conducted for 
maximum Reynolds number, Re max , of 50 and a Valensi number, Rc ()) , of 0.229. More details of this study 
are given in section 6.0. 

The results of this CFD research have been validated by comparing the CFD data with the literature 
and recent experimental correlations obtained at UMN (LSMU, large scale, friction- factor equation 3.5) 
and at Sunpower (actual-scale geometry, friction- factor equation 3.2 and Nusselt number equation 3.4). 
The agreement between these experimental correlations and the CFD data was good. 

For the steady-state 3-D simulation, both the local friction factor and the local-mean Nusselt numbers 
(i.e., mean value from the channel entrance to the local flow direction position) depart from the 2-D 
simulation values upon entering the second layer. That is where the 3-D effects become obvious and they 
persist as the axial coordinate advances. At the entrance of every layer, the forced reorientation of the 
flow results in small rises of both the friction factor and the mean Nusselt number with subsequent 
decrease as the flow settles into the new layer. Overall the plots of the friction factor and the mean 
Nusselt number tend to flatten out as the flow reaches a fully developed condition. 

As for the oscillatory flow simulations in 2-D and 3-D (with 6 layers), a base case was chosen using 
helium as a working fluid, stainless steel for the solid material, 3 10 K at the hot end, 293 K at the cold 
end, a maximum Reynolds number, Re max , of 50 and a Valensi number, Re co , of 0.229. Different 
parameters were examined (utilizing the 2-D model) to study the effects of changing: 1) the oscillation 
amplitude and frequency, 2) the Thermal Contact Resistance (TCR) between layers, and 3) the solid 
material. The effects of these parameters on the total regenerator heat loss (convection and conduction) 
were documented and are expected to be a useful tool for further development of Stirling engine 
regenerators. The baseline case with zero TCR, i.e., perfect thermal contact between layers, showed total 
losses of 2.896 W which was split between enthalpy losses (1.722 W) and conduction losses (1.174 W). 
As the TCR changed to infinity (perfect insulation between layers) the total losses decreased by 14%. 
This came about as a result of an increase in the enthalpy losses of 13.8% and a decrease in the 
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conduction losses by 54.7%. As we changed the solid material from stainless steel to nickel (while 
keeping an infinite TCR) we obtained an increase in the total loss by 3.8% due to a decrease in the 
enthalpy loss by 5.1% and an increase in the conduction loss by 36.3%. 

CFD simulations for the experimental jet penetration study done by UMN have been postponed to 
Phase 111. This study will include a slot jet and a round jet to examine the jet spread angle and penetration 
depth. 

3.6.5 Analytic Support (Gedeon Associates) 

3.6.5. 1 Adapting the Sunpower FTB for a Microfab Regenerator 

During Phase II the eventual goal of installing a microfabricated regenerator into the Sunpower FTB 
(frequency test bed) Stirling engine for testing was considered. Various options of increasing difficulty 
and cost required to adapt the engine to an involute-foil microfabricated regenerator with expected 
increasing levels of performance were identified. 

3.6.5.1.1 Background 

The Sunpower FTB Stirling convertor (engine plus linear alternator) sits in a test cell as a prototype 
for an Advanced Stirling Convertor (ASC) under development for NASA funded from a companion NRA 
program to the CSU microfab regenerator program. The FTB is a good choice for trying out a 
microfabricated regenerator because it is available and operates at a relatively low temperature (650 °C) 
compared to the ASC (850 °C). Parts are easier to make and operation does not require special high- 
temperature considerations. The goal is to adapt the FTB to use a microfabricated regenerator with 
minimal changes. 

During Phase I of the CSU-NRA contract some estimates were made of the performance advantage 
for a space power engine using a microfabricated regenerator (see app. C). During those studies there was 
the luxury of completely optimizing the engine to take advantage of the new regenerator. With the FTB, 
that freedom was not available. A mostly-fixed engine was available, allowing a swap-out of the random- 
fiber regenerator for a microfabricated regenerator and maybe a few other changes. 

Here are the things Sunpower was willing to change, listed in order of increasing difficulty and cost: 

1. Use either of two available heater heads, allowing two different regenerator lengths. 

2. Make a new heater head allowing for an increased regenerator length, with the same outer 
diameter and a new acceptor heat-exchanger insert with increased flow resistance (to increase 
pressure drop for displacer tuning reasons). 

3. Make a new heater head as above, but also with a different outer diameter to accommodate a 
“thinner” regenerator. 

Sunpower requested adherence to the following constraints: 

1. Keep the current displacer rod diameter so as to be able to use the existing piston/cylinder 
assembly for the new regenerator and thereby eliminate that experimental uncertainty. 

2. Restrict the regenerator length to no more than 60 mm to avoid displacer cantilever support 
problems. 

3.6.5.1.2 Sage 1-D Code Model 

This Sage code modeling was done prior to any actual testing of the involute-foil test regenerator, so 
the regenerator is modeled as a simple foil-type regenerator. The material was stainless steel and foil 
element thickness was fixed at 1 5 pm. The flow gap (between involute elements) was allowed to float as 
an optimized variable. The solid conduction empirical multiplier Kmult was set to 1 .0 as a conservative 
estimate pending a better understanding of the advantages of interrupting the solid conduction path. 

In order to maintain the displacer phase angle with a fixed rod, the overall flow resistance of the 
rejector + regenerator + acceptor must be maintained. So the lower flow resistance of the microfabricated 
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regenerator compared to the random fiber regenerator it replaces could not be used to full advantage. The 
model maintained flow resistance indirectly by imposing the constraint that the excess displacer drive 
power ( Wdis ) be zero. The displacer was constrained by a “displacer driver” component and Wdis is the 
power it requires to move the displacer at the desired amplitude and phase. To satisfy this constraint the 
model optimized the regenerator flow gap and, sometimes, the acceptor flow passage dimension. When 
the acceptor flow passage dimension was optimized there was enough slack to also maximize efficiency. 
Otherwise the regenerator flow gap was used only to meet the tuning constraint without regard to 
efficiency. 

3.6.5.1.3 Performance Estimates 

The simulated performances for the microfabricated regenerator with the above progressive FTB 
accommodations are given in table 3.5, compared to the baseline random fiber performance (via the ratio 
of the improved microfab efficiency to that of the baseline). 


TABLE 3.5— SIMULATED PERFORMANCE RESULTS 



W 

vv pv? 

w 

Qi„, 

w 

PV 

efficiency 

Efficiency/ 

baseline 

Baseline random-fiber regenerator 

128.0 

303.7 

0.422 

— 

Microfab — existing head, L reeen = short* 

123.8 

292.1 

0.424 

1.005 

Microfab — existing head, L reeen = long* 

113.2 

260.8 

0.434 

1.028 

Microfab — new head, same OD, L reBen = 60 mm, smaller acceptor passages 

101.7 

227.5 

0.447 

1.059 

Microfab — new head, smaller OD, L reg en = 60 mm, smaller acceptor 
passages 

110.7 

247.2 

0.448 

1.062 


* Short and long regenerators refer to two different dimensions allowed for design and not available for publication at this stage of the work. 


The final efficiency value falls in the range of the 6 to 9% efficiency gain projected earlier (see 
app. C). Most of the efficiency gain is achieved by making a new heater head, slightly longer but of the 
same diameter. The added bother of reducing heater head diameter hardly increases efficiency at all, 
though it does increase power level somewhat. 

The higher efficiency values of the last two micro fabrication cases (table 3.5) are a result of the 
acceptor taking the burden of providing the additional pumping dissipation to maintain the displacer 
phase angle, freeing the regenerator gap to optimize for efficiency. The acceptor is the best place to put 
the added dissipation because at the hot temperature it should, in theory, be partially recoverable. The 
extra dissipation subtracts from the available PV power, at-least somewhat, and it would be better to 
reduce the displacer rod diameter instead. But that is not allowed in the current FTB and will have to wait 
until the future when a space-power convertor might be designed from the ground up to employ a 
microfabricated regenerator. 

3.6.5.1.4 Regenerator Dimensions 

The regenerator flow gaps for the microfab cases of table 3.5 are given in table 3.6. 


TABLE 3.6— FLOW-GAP VALUES 



Flow gap, 
pm 

Microfab — existing head, L reeen = short* 

52.3 

Microfab — existing head, L reeen = long* 

58.0 

Microfab — new head, same OD, L reBen = 60 mm, smaller acceptor passages 

92.7 

Microfab — new head, smaller OD, L reBen = 60 mm, smaller acceptor passages 

91.6 


* Short and long regenerators refer to two different dimensions allowed for design and not available for publication at 
this stage of the work. 


The last two have much larger gaps (presumably much easier to make) as a result of the acceptor 
providing the additional damping to maintain the displacer phase angle. The overall regenerator 
dimensions are not included in the table to avoid revealing any proprietary FTB dimensions. 
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3.6.5.1.5 Recommendations 

Converting a FTB to use a microfabricated regenerator is feasible and it appears that the best option is 
the next to the last — a new longer heater head with the same OD and a re-designed acceptor insert. It 
remains to be worked out exactly how to hold the new regenerator in place and what measures will be 
required at either end to distribute the flow from the acceptor and rejector. 

3.6.5.2 Radiation-Loss Theoretical Analysis 

During the Phase 1 final review an evaluation of the effects of radiation through the regenerator was 
requested. This was done by a simplified theoretical analysis of radiation down a long thin tube (the 
results were later confirmed by a CSU CFD analysis). A long thin tube overestimates radiation through a 
stack of involute-foil disks because there is a clear sight path down the whole length of the tube. In the 
actual involute-foil stack it is impossible to see any light passing through it when held up to a bright light 
source. 

Looking down a long thin tube — eye focused to the far end — one sees mostly wall. So, too, radiation 
emitted at the hot end of such a tube sees mostly wall, where it is absorbed before it gets too far, provided 
those walls are diffuse gray absorbing surfaces (not highly reflecting). Emitted radiation gets about 3.5 
tube diameters before 99% of it hits the tube wall (easy to work out by comparing the surface area of a 
hemisphere with radius 3.5 tube diameters with that of the tube cross-section area). Some fraction of that 
radiation is reflected but it too cannot get very far before multiple reflections eventually absorb practically 
all of it. The walls of a long thin diffuse-gray tube act like a sort of distributed radiation shield. 

To the extent that the microfabricated regenerator (stack of involute-foil disks) looks like a bundle of 
long thin tubes it will also block radiation transmission. The analogy is not too far fetched. A view at the 
hot end of the regenerator looking toward the cold shows mostly foil surfaces, except for a tiny view 
angle where the cold end is visible. So, a quantitative estimate of radiation loss down a long thin tube can 
serve as a basis for estimating the radiation loss in a microfab regenerator. Perhaps it is not too 
unreasonable to substitute passage hydraulic diameter for tube diameter. 

The remainder of this section considers the limiting case of diffuse black regenerator walls, which is 
arguably a reasonable approximation for long narrow passages. The results are summarized in table 3.7 
which shows that radiation loss under these assumptions is negligible compared to other regenerator 
losses: 


TABLE 3.7— RELATIVE LOSS ESTIMATES FOR A 100 W 
CLASS SPACE-POWER STIRLING REGENERATOR 


Hot temperature 850 °C (1 120 K) 

Cold temperature 100 C (370 K) 

Passage aspect ratio L/d 300 

Passage wall emissivity s 0.5 

Radiation flow at cold end 10 mW 

Radiation flow at hot end 200 mW 

Time-average enthalpy flow 13,000 mW 

Solid conduction 7,000 mW 


3.6.5.2.1 Radiation in Long Thin Tube 

A tube of radius a (diameter d) and length L, as shown in figure 3.16, was used in this study. 
Coordinate 2, = x Id is the dimensionless axial coordinate. The tube has open ends at c, = 0 and tj/. = L/d. 
The tube wall is presumed to be a diffuse black surface (emissivity 8=1) with a fixed wall temperature 
T(g) that varies linearly with c between 71, and 7). The tube terminates in black cavities at the two ends at 
temperatures T 0 and 7). Of interest is the radiation heat flux q(q) through the tube section A(q). 
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To T l 



^o=0 \ q,= L/2a 


Figure 3.16. — Schematic of the tube used to study radiation heat transfer in the regenerator. 


Radiation Loss Multiplication Factor 
Long Tube, T 0 /T L = 1/3 



x/L 

Figure 3.17 — Radiation loss estimates. 



L/d = 10 
L/d = 30 
L/d = 100 
L/d = 300 


If radiation flux q(q) is small compared to the helium time-average enthalpy flux and solid thermal 
conduction flux down the regenerator (when installed in a running engine) then it will have a small effect 
on the usual regenerator temperature distribution and the assumption of a linear temperature distribution 
is valid. 

The radiation flux depends on position and can be represented as a fraction R of the worst-case 
radiation flux: 


= max (3-6) 

Where r/ max is the black-body radiation exchange between two parallel planes at temperatures 7'n and T L 
(flux limit as tube length approaches zero): 


q max =-a( T L- T o) (3-7) 

Constant a is the Stefan-Boltzmann constant 5.729E-8 W/(m 2 K 4 ). 

In general, multiplication factor R depends on the position q, temperature ratio T 0 /T L and tube aspect 
ratio L/d. Numerical calculation in a custom-written Delphi Pascal program for the representative case 
T 0 /T l = 1/3 (typical for Stirling engine) and various values of Lid gave the results shown in figure 3.17. 
L/df for the involute-foil regenerator is ~350 for a total regenerator length of 60 mm. In figure 3.17, the 
cold end is at x/L = 0.0 . 
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3.6.5.2.2 Justifying the Black-Wall Assumption 

Even though the actual wall emissivity for a metallic regenerator material is probably closer to s = 0.5 
than e=l, multiple reflections in a long thin tube render the apparent tube-wall emissivity near one at any 
location. The apparent emissivity is the value s a for which the total outgoing radiation (emitted + 
reflected) is e a o 7 4 (q). Figure 8.9 onp. 257 of (Siegel and Howell, 1981) shows that the apparent wall 
emissivity for an isothermal tube approaches 1 within a few diameters of the tube entrance, regardless of 
actual emissivity. In particular, for actual emissivity 8 = 0.5 the apparent emissivity is nearly 1 at a 
distance of only 2 tube diameters from the entrance. For the present analysis this means that it is not 
unreasonable to consider the tube walls to be black surfaces, which greatly simplifies the analysis because 
reflected radiation need not be considered. This black-wall assumption is arguably valid so long as the 
wall temperature does not change much over a distance of several diameters, which implies that the local 
radiation environment is similar to that of an isothermal tube. 

3.6.5.2.3 Applied to a Regenerator 

For the temperatures of a space-power engine, T L might be on the order of 1 120 K and T 0 might be on 
the order of 370 K (77/3) so the worst-case parallel-plate radiation heat flux q mm works out to 8.9 W/cm 2 . 
A regenerator void frontal area (corresponding to the tube interior of the above analysis) for a 1 00 W 
class Stirling engine is on the order of 2 cm 2 (see app. C). So the total worst-case radiation flow would be 
about 18 W. The actual radiation flow down the regenerator would be reduced by a fraction correspond to 
the curve R(q) in the above plot for Lid = 300. The aspect ratio for our current microfab regenerator 
design (based on length to hydraulic-diameter ratio) is about 350 for a total regenerator length of 60 mm. 

The conclusion is that the radiation heat flow down the regenerator would be very small. Near the 
cold end about 6x10 4 of 18 W, or about 10 mW. Near the hot end about 1x10 2 of 18 W, or about 
200 mW. In terms of the time average enthalpy flux (13 W) and the solid thermal conduction (7 W) the 
radiation loss is smaller by two orders of magnitude. For detailed derivation of the radiation heat transfer 
loss see appendix A. 

3.6.5.3 Solid Thermal Conduction in Segmented Foil Regenerators 

One of the selling points for our stacked-disk design was that it interrupts the solid thermal 
conduction path from one end to the other. This section explores this statement in detail and finds that the 
truth is more complicated and depends on disk thickness, solid thermal conductivity and also the 
properties and Reynolds number of the gas flowing through it. 

An analysis of segmented foil regenerators shows a complicated reality. In one extreme, coupling 
between the regenerator gas and solid bridges the contact resistance between segments, producing solid 
conduction in individual segments approaching that of a continuous foil regenerator. In the other extreme, 
high thermal conduction within each segment produces a stair-step solid temperature distribution with 
distinct temperature gaps between segments, increasing the net enthalpy flow down the regenerator. In 
either case solid conduction shows up as a regenerator thermal loss. The Mezzo regenerator lies closer to 
the second extreme, which is likely the cause of the reduced figure of merit at low Reynolds numbers in 
recent testing. This raises concerns about using high-conductivity materials, such as the pure metals 
(Nickel, Gold and Platinum) favored by the current electroplating fabrication process (FiGA). If such 
materials are used, the optimal regenerator disks need to be shorter and have thinner walls than disks 
made of lower conductivity materials 

3.6.5.3.1 Average Solid Conduction 

Equation (3.8) is a simple approximation for the average solid conduction in a segmented foil 
regenerator (Symbol definitions for this section are shown in table 3.8) 
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(3.8) 


lh_~ 

0,0 1 + 2 M )^1 

p k F 


where F is the not-quite so simple factor: 


2Nu L c h 2 

l-e ( RePr )A> (3.9) 

v ) 


TABLE 3.8— SYMBOL DEFINITIONS FOR THIS SECTION (3.6.53) 


Qs 

Spatial average solid conduction for segmented foil regenerator 

QsO 

Conduction for continuous foil regenerator with same solid area 

p 

Regenerator porosity (void fraction) 

A, 

Hydraulic diameter (twice flow gap) 

L c 

Foil segment length (disk thickness of Mezzo design) 

k s 

Solid thermal conductivity 

K 

Gas thermal conductivity 

< Rc Pr > 

Time averaged Re Pr product 

Rc 

Reynolds number based on hydraulic diameter 

Pr 

Prandtl number 

Nu 

Nusselt number ( hD h / k ) 


F = ( Re Pr)^- 
D h 


The physics behind this approximation is that heat-transfer coupling between the gas and solid tends 
to short-circuit the contact resistance between solid segments, reducing the temperature gap between 
segment endpoints. This imposes temperature gradients within the segments and therefore thermal 
conduction. The thermal conduction varies along the segment length but the average value is of interest. 
The details are derived in the following subsections of this section 3. 6. 5. 3. 

The point in looking at average solid conduction is that when it is high — approaching that of a 
continuous-foil matrix — then any contact resistance between segments is doing no good. Is that the case 
for the Mezzo regenerator? Table 3.9 shows some key values for that regenerator: 

TABLE 3.9— KEY VALUES FOR MEZZO 


INVOLUTE-FOIL REGENERATOR 


P 

0.84 

LJD h 

250/2(85) = 1.5 

kjk 

86/0.18 = 480 


For these values the average solid conduction ratio of equation (3.8) reduces to equation (3.10). 


Qs i 

Qs o ~ i + M 
F 


(3.10) 


The factor F depends mainly on Reynolds number and relative segment length LJD h and is plotted in 
figure 3. 1 8 for the case of Prandtl number, Pr = 0.7 and Nusselt number, Nu = 8.23 (developed flow 
between parallel plates, constant heat flux boundary condition). 

As relative segment length L c /D h increases, F grows quickly so that solid conduction approaches that 
of a continuous-foil regenerator. For the Mezzo regenerator with L c /D h = 1.5, the F factor peaks at about 
15, at a Reynolds number of about 30, which means that solid conduction is never more than about 8% of 
continuous foil conduction, according to equation (3.9). In other words the Mezzo regenerator segments 
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Solid Conduction Reynolds Number Dependence 



Figure 3.18. — Factor F of equation (3.9) as a function 
of Reynolds number at various segment lengths 
Lo/Dh, for Pr = 0.7 and Nu = 8.23. 

never have more than about 8% the regenerator-average temperature gradient. Each segment is relatively 
isothermal. The high contact resistance does indeed block most of the solid thermal conduction. 

3.6.5.3.2 Price of Localized Solid Conduction 

But isothermal regenerator segments are not without their cost. A stair-step solid temperature profile 
(piecewise constant) results in a minimum enthalpy flow per unit flow area of 

hmm=-c p p(u)KT g (3.11) 

where <u> is the time-averaged flow velocity (section-mean, absolute value) and A T g is the solid 
temperature difference between successive regenerator segments (see fig. 3.19). This is easy to 
understand by considering a regenerator cross-section between segments (disks) where the overall 
regenerator temperature is increasing in the positive direction. For positive flow the gas temperature is 
always lower than the solid temperature of the negative segment because it is heating up overall. For 
negative flow the gas temperature is always higher than the solid temperature of the positive segment. So 
the time-average difference in gas temperatures passing through the gap must be at least A T g . The 
minimum enthalpy flow can be put in dimensionless form by dividing by molecular gas conduction -k 
dT/dx. A key observation is that the overall regenerator temperature gradient dT/dx is just AT g /L c for a 
stair-step temperature distribution. After a bit of simplification the minimum enthalpy flow for a stair-step 
temperature distribution reduces to 


hmin ~ ~k ~~~ ^Rc Pr) ~~~ 
ax Dfr 


(3.12) 


Compare this to the following continuous-regenerator enthalpy flow, of equation (3.12), p. 52 of the 1996 
regenerator test-rig contractors report (Gedeon and Wood, 1996): 


dx \ 4Nu j 


(3.13) 


At low Reynolds numbers the minimum enthalpy flow h m in dominates h while the reverse is the case 
at high Reynolds numbers. The Reynolds number at which the two are the same is found by equating h m in 
to h and solving. That critical Reynolds number, assuming constant Nusselt and Prandtl numbers, is 
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(3.14) 


< Re l 


4Nu L c 
Pr D h 


Whenever Reynolds number is higher than this, the effect of the stair-step temperature distribution 
will be of diminishing significance. For the Mezzo regenerator <Re> c is about 50, where the value of h mm 
from equation (3.12) is about 50 times higher than the helium static conduction loss. If the regenerator 
were continuous foil the solid conduction loss would be a factor k s /k (l~P)/p higher than helium static 
conduction (conductivity ratio times area ratio) or about 91 times higher. So the segment-localized solid 
conduction produces an adverse helium enthalpy flow that is nearly as bad as if the solid conduction 
passed unimpeded by contact resistance through the entire regenerator. 

This coupling of localized solid conduction to helium enthalpy flow likely explains why the 
experimental test results for the Mezzo regenerator showed a suspiciously low Nusselt number below 
Re = 100 (as in fig. 3.14). The derived Nusselt number was forced to account for the h mm enthalpy flow 
based on a regenerator model that did not assume a stair-step temperature distribution. High solid 
conduction was at the root of the problem, in spite of the high contact resistance between regenerator 
disks. Had the solid conduction been lower, then the temperature gaps would have been smaller and h mm 
would have been smaller. 

What might be done about this problem? The obvious solution is to make the regenerator disks from 
lower conductivity material. The ratio of solid to gas conductivity k s /k is rather high for a nickel 
regenerator, or about 5 times that for a stainless-steel regenerator. As a rule of thumb, alloys have lower 
thermal conductivity than pure metals. The following table (table 3.10) gives some idea of the solid 
conduction losses of pure-metal regenerators compared to a stainless-steel regenerator. 


TABLE 3.10— SOLID CONDUCTION 


Material 

Thermal conductivity near room 
temperature* (W/m C) 

Ratio relative to 316 
stainless steel 

316 stainless steel 

16 

1 

platinum 

73 

4.6 

nickel 

86 

5.4 

gold 

300 

19 


* Source: 1985 Material Selector Handbook, Materials Engineering 


Let’s not be too pessimistic here. The good figure-of-merit measured for the nickel regenerator 
includes the effects of solid conduction loss. But the involute-foil would do even better in the future with 
a lower conductivity material or thinner walls, especially at lower Reynolds numbers where the 
temperature-gap effect dominates. The Mezzo regenerator figure of merit peaked at a Reynolds number of 
about 200 whereas most Stirling engine optimizations performed in the past prefer Reynolds number 
amplitudes on the order of 1 00 or lower. 

3.6.5.3.3 Solid Conduction Derivation 

The solid temperature in a segmented-foil regenerator lies somewhere between the two extremes 
illustrated in figure 3.19 — between a piecewise constant stair-step distribution and a continuous 
distribution. The first case corresponds to static-conduction in vacuum with very high contact resistance 
between segments. There is a temperature gap A T g between segments. In the second case the gas is in 
very good thermal contact with the solid and has sufficient heat capacity to dump any required amount of 
heat in the segment entry regions so as to eliminate the temperature discontinuities. 

The reality is somewhere between the two. That reality looks something like the illustrations in 
figure 3.20 which shows the gas and solid temperature distributions (time averages) for the middle two 
segments of a 6 segment (shorter) segmented regenerator. The solutions were generated by a 1 -D 
Sage-code simulation, as discussed in detail later. Gas temperatures are averaged over the positive and 
negative flow half-cycles individually. Without time averaging, the temporal temperature variations (due 
to finite solid heat capacity) would appear significant at this scale, but they occur roughly 90° out of 
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Figure 3.19. — Extreme limiting cases for solid temperature distribution in a 
segmented foil regenerator. 


phase with the mass flow rate; so they do not affect the net enthalpy flow. The upper illustration shows a 
solution for a segmented regenerator with the same properties as the involute-foil regenerator (250 pm 
thick), except with stainless-steel solid material in order to better show the solid temperature variation. 

The lower illustration is for a stainless-steel regenerator with 3x longer segments, corresponding to 
750 pm thick disks, and a higher Reynolds number, where the solid temperature gradient increases to 
about 60% of the regenerator average. 

In any case the average solid conduction within a segment is proportional to the average solid 
temperature gradient, which may be written (A T m - A T g )/L c , where AT,,, is the temperature difference 
between neighboring segment centers, A T g is the temperature gap between neighboring segment 
endpoints and L c is the segment length. The average solid conduction is then the product of solid 
conductivity k s , solid cross-section area A s and the average temperature gradient: 

Qs - ~k s A s — (3.15) 

L c 

The solid temperature gap A T g depends on the overall energy balance which can be roughly 
formulated in terms of a control volume between the beginning of a regenerator segment (left) and its 
middle (right), as shown in figure 3.21. 

The energy balance is formulated in terms of the time-average enthalpy flow H carried by the gas and 
the axial solid conduction Q s . At the left boundary (segment endpoint) these values are H e and zero, 
respectively (zero solid conduction because of presumed high contact resistance between segments). At 
the right boundary (mid segment) the values are H, and Q sc . The half-segment energy balance may 
therefore be written as in the following equation (3.16): 

H e= H c+Qsc (3.16) 

The heat transfer q between gas and solid, shown in figure 3.21, is not part of the energy balance because 
both gas and solid are included together and what leaves one enters the other. 

The essential idea is that the inter-segment temperature gap A T g results in an increased net enthalpy 
flow at the segment entrance compared to the equilibrium value deep within the segment after the gas has 
had time to drop its extra heat to the solid. For a long segment the increased enthalpy flow is —cp< m > 

A T g , where < m > is the time-averaged absolute mass flow rate. For a short segment the increase is not as 
much. In a section of any length the difference in net gas enthalpy flow between the segment end and its 
midpoint is approximately: 
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(3.17) 


H e -H c 


-c, 
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where the extra factor (1 - e A ’) is derived later under the sub-title, “Steady Heat Transfer Solution,” with 
A given by equation (3.23). 


Sage Solution Segmented Foil Regenerator 



Sage Solution Segmented Foil Regenerator 
3x Longer Segments 



Figure 3.20. — Time-average solid- and gas-temperature distributions 
for the middle-two segments of a 6-segment foil regenerator. Top: 
250 pm long (shorter) segments, Bottom: 750 pm long (longer) 
segments. 
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Figure 3.21 . — Energy balance between segment end (left) and center 
(right), bounded by vertical dotted lines. The solid curve represents 
time-average solid temperature. The dotted curves above and 
below represent gas temperatures time-averaged separately for 
positive and negative flow directions. 


The dumping of heat from the gas to the solid results in solid conduction that varies smoothly from 
zero at the segment endpoint to a maximum at the midpoint. If the segment is not too long then it is 
reasonable to assume that the average solid conduction is halfway between the endpoint and midpoint 
values. In other words: 

Qsc=2Qs (3.1 

In terms of the above simplifications the energy balance can be written as in equation (3.19): 


-c p (th)AT g l-< 


e /2 *2Q S 


Solving for A T g and substituting into equation (3.15) gives the average solid heat flux in the form 


— AT k A 

f) ~ k A m ? s s 
\£s ~ As T A T 


, (iii) 1 - 


The first term on the right is just the conduction Q s0 for a continuous foil with the same solid cross section 
and boundary temperatures. The second term contains several factors that can be arranged into more 
standard dimensionless groups. Solving for Q s the result is equation (3.8), given near the beginning of 
this section. 

The average segmented regenerator conduction approaches the continuous foil conduction as the term 
2 — — — ^ of equation (3.8) approaches zero, which can happen for a number of reasons. For 

example, high porosity (3 (thin foil), low solid conductivity k s , long segment length L c , or low Reynolds 
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number. On the other hand the average segmented regenerator conduction approaches zero as that same 
term approaches infinity, which happens for the opposite reasons. 

3.6.5.3.4 Sage Models 

The illustrations in figure 3.20 were the result of modeling individual segments of a foil regenerator 
using the Sage code. It would be very awkward and time consuming to model a total regenerator this way 
(consisting of hundreds of segments) but it is not too difficult to model, say, a 6-segment regenerator in 
order to get an idea of what is going on. The Sage model illustrated in the following figure 3.22 does just 
that. 

This Sage model is equivalent to the NASA/Sunpower regenerator test-rig containing a rather short, 
6-disk, involute-foil regenerator. Each canister of the model contains a simulation of a single involute-foil 
regenerator disk, identical to one fabricated by Mezzo except made of stainless steel instead of nickel. 

The reason for stainless steel is because its lower thermal conductivity results in a higher solid 
temperature variation which shows up better on plots. The regenerator segments are inter-connected by 
gas-flows (including gas thermal conduction continuity) but the solid domains are not connected, 
corresponding to high contact resistance between disks. 

Middle segments A and B, in figure 3.22, are the segments of main interest (whose solutions are 
plotted in fig. 3.20) with the “buffer” segments on each side included to give some opportunity for the 
temperature solutions to develop to their equilibrium (spatially periodic) values. 

3.6.5.3.5 Baseline Sage Model 

The working gas is helium at 25 bar charge pressure and the driving piston motion is sinusoidal with 
amplitude adjusted to produce an average Reynolds number of 30 in the regenerator segments. This 
Reynolds number gives the largest F factor according to figure 3.18. The endpoint boundary temperatures 
are adjusted so the temperature difference between segments is about 2.5° corresponding to about what it 
would be in a realistic Stirling-engine regenerator (e.g., 800 °C temperature difference over 320 disks). 
The segment temperatures range from about 293 to 308 K, so the model corresponds to a short section 
near the cold end of the regenerator. 



Figure 3.22. — Sage model diagram for 6-segment regenerator. 
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The Nusselt number in each segment is set to the constant value 8.23, corresponding to developed 
laminar flow between parallel plates (constant heat-flux boundary condition). This is the most unrealistic 
part of the model since it neglects the increased Nusselt number in the developing-flow region, which is 
arguably the entire section. But the intent of the model is only to get a rough idea of what is going on for 
purposes of validating the simplified theory derived in this memo. It does not make sense to use the 
Nusselt number derived from recent testing of the Mezzo regenerator because that Nusselt number applies 
to the regenerator as a whole and not to individual segments of the detailed model. 

3.6.5.3.6 Longer Segments 

The baseline model is not very satisfying from an academic viewpoint because there is not much solid 
temperature variation within segments. This can be remedied by increasing the segment length by a factor 
of 3, producing a ratio L c /D h = 4.5. To keep the same overall temperature gradient the temperature 
difference between segments is also increased by a factor of 3. By increasing the average Reynolds 
number to 100 (increasing charge pressure and stroke) the F factor of figure 3.18 is roughly at its 
maximum value of about 100 and the average temperature of the solid should be about 70% of the 
regenerator average according to equation (3.8). That is not too far off from the 60% of the Sage solution 
shown in the lower illustration of figure 3.20. 

3.6.5.3.7 Steady Heat Transfer Solution 

At the heart of the preceding solid conduction derivation is the energy balance in a single segment of 
the regenerator. That energy balance is the result of a simplified steady-state solution for heat transfer 
between a solid segment of length L c with a time- invariant linear temperature distribution T,,(x) = mx and 
two steady gas streams of mass-velocities ±pw, as illustrated in figure 3.23. 

The governing equation, equation (3.21), is a simplified energy equation for incompressible, steady 
flow, considering only heat transfer to the solid: 

~~ — ~~ {T ~T X ) (3.21) 

ax Cppu 


In figure 3.23, the gas temperature solution for positive directed flow is denoted T+ and the solution for 
negative flow 71. The boundary condition is that there is a temperature difference A 71 between the gas 
and solid at the negative end for T+ and the positive end for 71. 

Of interest is the temperature difference 71 - 71 at either entrance, denoted A 71, compared to the 
temperature difference at the segment midpoint, denoted A T c . That temperature difference determines the 
amount of heat transfer to or from the solid between the segment end and midpoint. Skipping all the 
details, the final result is: 
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where A is: 
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Steady Temperature solutions in Regenerator 
Segment 



T 

A T 0 


relative position in segment 

Figure 3.23. — Steady-flow gas-temperature solutions for regenerator solid 
segment with linear temperature distribution. 


The asymptotic temperature difference T-T s for a long section is just ±mL c IA, depending on the flow 
direction. In that case the first factor on the right in equation (3.22) is the total gas temperature change 
from section entrance to exit, which is the same as the solid temperature difference A T g between 
regenerator segments. Making this approximation regardless of segment length results in the following 
approximation, equation (3.24), for the change in net gas enthalpy flux (per unit flow area) between the 
ends and middle of a regenerator segment: 

h e - h c &-c p p\u\/iTg{\- e~ /l\ (3.24) 

This is essentially the approximation used in energy balance equation (3.17). In making the leap from 
steady flow to sinusoidal flow the factor \u\ appears in equation (3.17) in the form of the time-average 
absolute mass flow rate, the time variation of the factor ( l—e A 2 ) is ignored and A is evaluated at the mean 
Reynolds number. 

3.6.6 Structural Analysis of Microfabricated Involute-Foil Regenerators (Infinia) 

Finite element analysis of the microfabricated involute-foil regenerator showed that the regenerator 
has very high average axial direction stiffness 6.5xl0 9 N/m (3.75 x ] Q 7 lb/in.). Without any radial side 
disturbance, the stress level was much lower than the material yielding strength. If a radial side 
disturbance such as misalignment is localized in a small area, the Von Mises stress is beyond the material 
yielding strength and permanent deformation could occur in that area, which may decrease the Stirling 
efficiency. In order to prevent local permanent deformation, the radial side load must be small, or the 
disturbance area must be large. 

In summary, the proposed microfabricated involute-foil regenerator has high axial stiffness. The 
stress level is sensitive to the radial side disturbance, which therefore requires special cautions and an 
appropriate process for installation to prevent lateral permanent deformation. 

The structural analysis is more thoroughly discussed in section 7.0. 
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3.7 Summary of Expected Benefits of Segmented-Involute-Foil Regenerators 

The power conversion efficiency obtained with this new regenerator is estimated to show a 6 to 9% 
overall improvement based on calculations for the Sunpower ASC engine. The specific power (We/kg) 
will be directly proportional to the efficiency. As for lifetime and reliability, the microfabricated 
regenerator has the potential to remove any concerns about stray fibers, as may result from use of 
random-fiber or wire-screen matrices. An estimate of development risk is “medium” since the 
manufacturing viability has been established but the concept has not yet been tested in an engine. 
However, an involute-foil-layered regenerator has been very successfully tested in the NASA/Sunpower 
oscillating-flow test rig (with results shown in fig. 3.9). And involute-foil regenerator microfabrication is 
currently underway for testing in an engine during the Phase III effort. 

4.0 Mezzo Microfabrication Process 

4.1 Background 

This section of the report focuses on the contributions of International Mezzo Technologies. One 
important reason Mezzo was selected as the subcontractor for this project was its expertise with the LiGA 
micromachining process and its investment and commitment towards combining the advantages of LiGA 
with Electric Discharge Machining (EDM). The combination of LiGA-EDM theoretically provides a 
means to fabricate high aspect ratio micro features normally associated with the LiGA process, out of any 
conducting material. 

The initial research plan involved Mezzo using the LiGA process to fabricate well-defined, high- 
aspect-ratio EDM tools. These LiGA-fabricated EDM tools would then be used to make the 
micromachined regenerator parts from materials with the desired high temperature properties and low 
thermal conductivity (stainless steel, Inconel, etc.). EDM tools were fabricated via LiGA and efforts to 
EDM parts from stainless steel showed initial promise in terms of being able to produce the correct 
geometry, at least at shallow depths. But the process was very slow, tool wear rate was high, and it 
became apparent that the probability of fabricating the desired stainless-steel regenerator using LiGA- 
EDM with the available funding was low. 

To fabricate the regenerator on schedule, Mezzo changed its manufacturing approach. The standard 
LiGA process was used to directly produce individual nickel regenerator components which were then 
assembled, sent to Cleveland State University, and subsequently tested at Sunpower. By changing its 
fabrication strategy, Mezzo was able to provide the regenerator for the project and Sunpower was able to 
experimentally verify that the involute-foil-layered regenerator geometry provided high performance. 

This change in plans also supported the desire to move the oscillating-flow rig testing from Phase III 
(year 3) to Phase II (year 2) — since available Phase III funding for this and other NRA contracts was in 
jeopardy. This section provides a summary of the successful effort that resulted in the manufacture via the 
LiGA process of the electroplated nickel regenerator that was tested by Sunpower in the oscillating-flow 
test rig. 


4.2 LiGA-EDM 

The need to change from a LiGA-EDM approach to a pure LiGA approach was initially supported by 
test results at Mezzo that focused on limiting overbum to what was believed to be an acceptably low 
10 pm. The ramifications of this assumption will be revisited in this section of this report. Mezzo used a 
Mitsubishi Model EA8 die-sinker EDM machine to define features. It was found that the combination 
providing the lowest power and highest frequency (setting El 01 3) provided the best definition of micro 
features with a low overbum (approximately 1 0 pm of overbum). At this setting, a feature on the EDM 
LiGA tool with width 60 pm would produce a corresponding feature width of 80 pm in the material being 
machined. Using that setting, material removal rates have been quantified. Table 4.1 provides the material 
removal rate of a flat LiGA tool having electrode area 72.77 mm 2 . The material removal rate of a variety 
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of materials was quantified. In these experiments, no piezoelectric stack was used to provide enhanced 
flushing. For reference it has been found that the use of the piezoelectric stack increased the allowable 
bum rate of erbium by a factor of three, so similar improvement might be expected with other workpiece 
material selections. Electrode wear, bum depth, volume removed, bum time, and removal rate were all 
measured and are provided in table 4.1. 


TABLE 4.1— EDM RESULTS USING LiGA TOOL (AREA = 72.8 mm 2 , PIEZOELECTRIC STACK OFF) 


Electrode 

material 

Workpiece 

material 

Machine 
bum setting 

Electrode 

wear, 

|im 

Burn 

depth, 

pm 

Volume 

removed, 

mm 3 

Bum time, 
hr 

Removal 

rate, 

mm 3 /hr 

LiGA Ni 

420 SS 

E1013 

15 

12.7 

0.92 

18 

0.05 

LiGA Ni 

Titanium 

E1013 

38 

38 

2.77 

24 

0.115 

LiGA Ni 

Aluminum 

E1013 

5 

33 

2.4 

18 

0.133 


It can be seen that the removal rate of 420 SS (stainless steel) is 0.05 mm 3 /hr. To put this removal 
rate in perspective, the nickel involute-foil regenerator consisted of approximately 45 disks, each around 
25.4 mm in diameter, each 0.25 mm thick. The total volume of the regenerator would then be 5700 mm 3 . 
Of that volume, if the assumption is made that the overall porosity is 80%, then the volume of material to 
be removed is 4560 mm 3 . To fabricate a similar regenerator from stainless steel, based on the material 
removal rate provided in table 4.1, would require 91,207 hr (10.4 years), not counting the time to change 
tools, etc. It is reasonable to assume that the use of the piezoelectric stack might reduce this time by a 
factor of two. Also, Mezzo does not use a micro EDM which allows combinations of higher frequency 
and lower power than does the more traditional die sinker EDM used at Mezzo. However, conversations 
with users of micro EDM equipment indicated that the material removal rates achieved at Mezzo were 
roughly comparable with typical micro EDM results. In summary, no reasonable scenario currently exists 
whereby, if the overburn is to be limited to a value on the order of 10 /urn, the time to produce the CSU 
regenerator would be less than 1 to 2 years. 

For reference, it is also worthwhile to compare the material removal rates achieved in this project 
with the potential material removal rates that can be achieved by the Mitsubishi machine at higher power. 
The maximum removal rate of stainless steel is approximately 2.5 in 3 /hr. This corresponds to a removal 
rate of approximately 4 1 ,000 mm 3 /hr. The material removal rate at the highest setting therefore exceeds 
by a factor greater than 800,000 the material removal rates associated with the Cleveland State project. 
However, the overbum associated with such high material removal rates is about 400 pm — far in excess 
of the allowable limit for the CSU involute-foil regenerator. 

The reasonable question needs to be asked: Is there an intermediate acceptable balance between 
overbum, LiGA feature width, and material removal rate. Recently, a preliminary set of tests was 
performed where the power settings were increased and the material removal rate was measured as a 
function of overbum. The goal of the tests was to plot for a given tool geometry the overbum versus the 
material removal rate. The result of one test is shown in table 4.2. The LiGA tool consisted of a 500 pm- 
wide LiGA nickel rib. Material removal rates and overbum are shown as a function of power setting. 


TABLE 4.2— MATERIAL REMOVAL RATES/ OVERBURN FOR 500 pm- WIDE LiGA FEATURE 


Power 

setting 

Electrode 
surface area, 
mm 2 

Workpiece 

material 

Volume 

removed, 

mm 3 

Bum 

time, 

hr 

Removal 

rate, 

mm 3 /hr 

Over- 

bum, 

pm 

Removal 
rate ratio 

CSU regen. 
burn time, 
yr 

E1013 

72.77 

420 SS 

0.9242 

18 

0.05 

9 

1 

10 

E1014 

2.806 

420 SS 

0.2707 

1.9 

0.14 

11.5 

2.8 

3.6 

E1015 

2.954 

420 SS 

0.3013 

0.2 

1.5 

18 

29.4 

0.34 

E1018 

2.9654 

420 SS 

0.6026 

0.09 

6.7 

32.5 

130.5 

0.07 

E1019 

2.9654 

420 SS 

0.5272 

0.08 

6.5 


127 

0.08 
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c) E1018 setting: channel width (5x) d) El 1 8 setting: bottom of channel (lOx) 

Figure 4.1 . — Comparison between power setting and surface finish of bottom (sides not seen). 


Comparing settings E1013 and E1018, it is possible to deduce that if the acceptable overburn is 
increased from 9 to 32 pm, the removal rate increases by 100 fold. In retrospect, Mezzo probably 
concentrated excessively on extremely low power settings that were mandated by another LiGA-EDM 
project. In that project, LiGA nickel posts of diameter approximately 30 pm were used to produce 
features (holes) of diameter 50 pm. No more than 10 pm of overbum could be tolerated. The same 
assumption seems to have been made in accessing the potential of using LiGA-EDM to fabricate the CSU 
regenerator. In retrospect, a tool width of 40 pm, at high power setting, might produce a 1 00 pm-wide 
channel at a material removal rate that would be acceptable. Figure 4. 1 shows EDM results for some of 
the EDM settings shown above. E1018 has a material removal rate that is 130.5/2.8 = 46.6 times faster 
than E1014. The recast region on the sidewall does not look bad. The RMS roughness of the bottom of 
the channel is yet unquantified, but higher for the E1018 setting 

The conclusion of the recent EDM effort is that if extremely small overbum is necessitated, then 
EDM may be unacceptably slow. However, if the relationship between EDM material removal rate, 
overbum, and acceptable regenerator geometries is known, there may be a design space where the 
combination of EDM and LiGA is a useful manufacturing option that should be explored. 


4.3 LiGA-Fabricated Regenerator 
4.3.1 Manufacturing Process 

Two closely related LiGA processes are described below. One is referred to as the “optimal” process 
that was originally envisioned, the other describes the process that was actually followed. The difference 
between the two processes is associated with unanticipated problems in the development component of 
the lithographic patterning of the polymethyl methacrylatete (PMMA, or Plexiglass type plastic) 
templates. 

In an optimum LiGA exposure-development sequence, the sidewalls of the lithographically- patterned 
PMMA template are straight, as shown in figure 4.2a. In this project, unexpected difficulties developing 
the PMMA resulted in excess material removal, or “undercutting,” at the PMMA-substrate interface as 
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shown in figure 4.2b. This undercutting leads to a more complicated two-step electroplating process and 
an extra step involving EDM. Because of “undercutting” during the development process, the following 
process (described in fig. 4.3) was used to fabricate the regenerator components: 

1 . An x-ray mask was fabricated. The mask consisted of a tightly packed array of nineteen 
regenerator disks patterns. 

2. A 250 pm-thick sheet of PMMA was bonded to a 400-series (magnetic) stainless steel substrate. 

3. An x-ray lithography-electroplating process sequence was used to produce the nickel regenerator 
disk parts. It was found that development of the exposed PMMA caused some unexpected, 
undesired “undercutting” at the PMMA-substrate interface. Undercutting is associated with 
excess PMMA being dissolved during the development process. This fact motivated a two-part 
electroforming process. A copper electrodeposition step was used to fill the bottom of the features 
with copper to a depth equal to the height of the “undercutting” region. Beyond this point, nickel 
was deposited. To ensure that all these voids were completely filled with metal, the electroplating 
process was continued after all the features were filled, resulting in an “overplated” deposit. 

4. Initially, polishing was tried to remove the overplated layer. The polishing was found to destroy 
parts, so an alternative process was used and found to be successful. This successful process 
involved attaching the conductive substrate to a magnetic chuck, orienting the substrate in the 
vertical plane (fig. 4.3a and b). Then a wire EDM was used to take a “skim pass” just above the 
non-conductive PMMA layer (fig. 4.3b). This step removed the overplated nickel. The substrate 
was then released from the chuck. 

5. At this point, the nickel and copper electrodeposited features and the remaining PMMA was 
debonded from the substrate (fig. 4.3c) and the unexposed PMMA was dissolved in acetone 
(fig. 4.3d). The remaining nickel-copper features were again attached to the magnetic chuck with 
the nickel features in contact with the chuck (fig. 4.3e and f). A second EDM process was used to 
remove the copper and nickel in the “undercut” region. It should be noted that if the copper had 
completely filled the “undercut” region, the copper could have been removed with an etch, 
leaving only nickel parts with the desired geometry. However, it was found that insufficient 
copper was deposited to fill the “undercut” region. As a result, some nickel was also deposited 
into the “undercut” region, making it necessary to use a second EDM “skim cut” to remove both 
the nickel and copper within the “undercut” region. Following the second EDM “skim cut” the 
parts were released from the chuck (fig. 4.3g) and inspected. 


PMMA 



“Undercut” region 


substrate 


b) PMMA when undesirable “undercutting” occurs 
during the development process 


Figure 4.2 — Comparison of desired and actual lithographic patterning of PMMA. 
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a) Plated sample as removed from the plating bath. b) Same sample mounted on a magnetic chuck in wire 
The polymer resist is still intact. EDM machine. The wire EDM machine skims any over 

plated nickel to the same surface as the polymer resist. 



e) Sample mounted to magnetic chuck with copper f) Second EDM pass removes copper and nickel 
plating exposed. associated with “undercut” region. 



g) Part separated from chuck, ready for inspection. 

Figure 4.3 — Actual process used to fabricate regenerator disks. 
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4.3.2 Fabricated Regenerator 

The previously described process was used to fabricate the regenerator tested in this project. 
Micrographs of typical parts are shown in figure 4.4a to d. The nickel webs are approximately 15 pm in 
width, and arranged in an involute pattern (fig. 4.4a and b). The thickness of each disk is approximately 
250 pm. Figure 4.4c shows a single involute-foil slipped onto the stacking fixture. Figure 4.4d shows a 
single disk leaning against the outer housing of the regenerator. Figure 4.5 shows the final regenerator 
that was tested. 



c) Disk stacked onto fixture. d) Disk leaning against outer housing. 


Figure 4.4. — Different magnified views of regenerator disks. 



Figure 4.5. — Assembled regenerator (stack of 42 disks). 
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4.4 Microfabrication Process Conclusions 


In this project, the LiGA micromachining process was used to fabricate a regenerator that was tested 
and found to provide very good performance (see sec. 3. 6.2. 4 for test results). Also, the manufacturing 
approach of using LiGA-fabricated EDM tools to fabricate regenerator parts seemed initially to offer little 
potential due to the extremely low material removal rate. In retrospect, however, if dimensions are chosen 
such that greater overbum is allowable, then the EDM-LiGA approach will be more viable. It should also 
be noted that in this effort, a new manufacturing technique was developed: namely using EDM to “skim 
cut” regenerator parts. While the regenerator tested did provide good performance, LiGA and/or LiGA- 
EDM process optimization could result in a better product. Potential improvements include: 

i) Improve the lithography process to eliminate or greatly reduce “undercutting.” 

ii) Cease electroplating before overplating begins. This would eliminate the need to use the “skim 
cut” and would eliminate the source of burrs that was attributed to the EDM operation. 

iv) To greatly reduce cost, explore the use of SU-8, a negative resist that requires substantially less 
exposure time than PMMA. 

iv) For EDM, find acceptable combination of material removal rate, overbum, and geometry that 
gives a high quality part in a reasonable time and cost. 


5.0 Large Scale Mockup (University of Minnesota, UMN) 

5.1 Large Scale Mockup (LSMU) Design 

The microfabricated regenerator is of an annular design that cannot be scaled up in its entirety by a 
factor of 30 and still be operational in our oscillatory flow facility. Thus, only a 30° sector of it was 
chosen for modeling. Two geometries are shown in figure 5.1a and b. The second pattern geometry is 
achieved by shifting and flipping the first geometry. Figure 5.2 shows the dimensions of typical LSMU 
channels. The channel width is 2.58 mm and the fin thickness is 0.42 mm. The channel length changes 
from 54.56 to 32.72 mm as one passes from the inner radius to the outer radius of the first pattern layer. 
For the channels of the inner and outer edge of the second pattern layer, the lengths are 22.43 and 
14.32 mm, respectively. For the channels of the center area, the lengths change from 49.73 to 34.65 mm. 
The layers were stacked to make the LSMU assembly. Figure 5.3 shows the mesh of the first pattern layer 
stacked on top of the second pattern layer. 



Figure 5.1(a) — First pattern geometry (6 ribs), (b) Second pattern geometry (7 ribs). 
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The LSMU regenerator layers were fabricated by wire EDM. Figure 5.4 shows a photo of a LSMU 
layer, which is a 30° sector of the first pattern design. The choice of aluminum was influenced by the 
EDM manufacturing choice for more rapid cutting compared to stainless steel, another reasonable choice. 
By inspection, the surface appears smooth with a matte finish. Literature indicates the roughness obtained 
by wire EDM is 0.05 pm. 



Figure 5.2 — Geometry of typical channels. 



Figure 5.3. — The mesh of the first pattern layer stacked on top of 
the second pattern layer. 
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Figure 5.4. — The LSMU layer fabricated by wire EDM. 


Top dead center (TDC) trigger 



Figure 5.5. — Schematic of the oscillatory-flow generator (Seume et al. 1992). 


5.2 The Operating Conditions 

For engine operational data on which to base the LSMU design, a “pattern engine” was selected. It is 
representative of modern Stirling engines being developed for NASA space power applications. The 
pattern engine was designed for operation with a random wire regenerator. Operational data for this 
pattern engine and the hydraulic diameter of the microfabricated regenerator were used to calculate the 
Reynolds number and Valensi number for our pattern engine with a microfabricated regenerator installed. 
The computed Reynolds numbers varied from 19.7 to 75.7 for the hot end to the cold end of the 
regenerator while the Valensi number varied from 0.12 to 0.6. Then the stroke, piston diameter and 
frequency of the oscillatory-flow test facility were selected for use in the large-scale experiment that 
matched these dimensionless numbers. A Scotch-yoke mechanism is employed to produce precise 
sinusoidal movement of the piston with a zero-mean velocity in the oscillatory-flow test facility. 

Figure 5.5 shows a schematic of the oscillatory-flow generator, the details of which were given in a 
NASA report (Seume et al. 1992). Figure 5.6 shows a picture of the oscillatory-flow generator. 

Figure 5.1a shows the geometry of the first pattern regenerator layer (with 6 ribs), where the outer 
radius, R a , is 284.25 mm and the inner radius, R„ is 77.25 mm . The second pattern regenerator layer 
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Figure 5.6. — The picture of the oscillatory-flow generator. 


(7 ribs) has the same outer and inner radii as the first pattern. The scale factor is 30 times the actual size. 
For the first pattern regenerator layer, scaled-up channel lengths change from 54.56 mm for the slots near 
the inner radius to 32.72 mm for the slots near the outer radius. Scaled-up channel width is 2.58 mm 
(101 mils). The hydraulic diameter of the channel, d h , is 4.87 mm (192 mils). The wall thic kn ess is 
0.42 mm (17 mils). The plate thickness is 7.95 mm (312 mils or 5/16 in.). The porosity is 86%. The area 
of the 30° sector of an annulus is n x (Ro 2 - Ri 2 ) x 30/360 = 0.01959 m 2 . 

The stroke, piston diameter and operating frequency are selected to match the Valensi number 
and Reynolds number of the pattern engine. For air, at ambient pressure and temperature, viscosity, 
v = 15.9xl0' 6 m 2 /s. The stroke and piston diameter are set to the most suitable values of the options 
available in the present rig: 

Stroke = 178 mm (7 in.); piston diameter, D p = 216 mm (8.5 in.). 

Assuming incompressible fluid and balancing the volumetric flows in the regenerator and the driving 
piston/cylinder zone: 


226 max x A x (|) — 


7t Di 


- x Stroke 


(5.1) 


The amplitude of fluid displacement based on the LSMU regenerator flow area is 


3 % 


7iZ)2 X Stroke rt0 .216 2 x 0.178 


8 x J x (|) 8x0.01959x0.86 

Operational frequency is 0.2 FIz. Thus, 


= 0.1936 m (7.6 in.) 


(5.2) 


Re r 


U ma x x d/ 7 _ X max x(nxdh _ 0.1936 x2xjix0.2x 4.87 x 10 -3 
v v 15.9xl0- 6 


= 74.5 


(5.3) 
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(5.4) 


y a - d JL^. - d ^ x 3.14 x 0.2 x (4.87 x 10~ 3 ) 2 _ ^ 

4v 4v ~ 4x15. 9xl0- 6 

The microfabricated regenerator (in the operating engine) Reynolds number varies from 19.7 to 75.7 
from the hot end to the cold end of the regenerator, while the Valensi number varies from 0.12 to 0.6. In 
our large-scale experiment, the Reynolds number is 74.5 and the Valensi number is 0.47, a reasonable 
match to the dimensionless numbers for the pattern engine. The local maximum velocity, U max , is 
0.24 m/s. Adolfson (2003) found that quasi-steady velocity measurements with hot-wire anemometry 
could be made with less than 10% uncertainty when the velocity exceeds 0.12 m/sec. 

Thus, in the large scale experiment, the stroke is 178 mm, the piston diameter is 216 mm and the 
frequency is 0.2 Hz, selected to match the Reynolds number and the Valensi number of the microfab 
regenerator in the pattern engine. Table 5.1 shows a comparison of microfabricated regenerator and 
LSMU regenerator geometry and working conditions. 

5.3 The LSMU Experiments Under Unidirectional Flow 

The Darcy friction factor, the permeability and the inertial coefficient of the LSMU layers were 
measured under unidirectional flow. 


TABLE 5.1— COMPARISON OF MICROFABRICATED AND LSMU INVOLUTE-FOIL REGENERATORS 


| Microfabricated regenerator | LSMU regenerator 

Geometry 

Channel width (mm) 

0.086 

2.58 

Channel wall thickness (mm) 

0.014 

0.42 

Regenerator layer thickness (mm) 

0.25 

7.9 

Hydraulic diameter, d/, (mm) 

0.162 

4.87 

! Working conditions j 

Working medium 

Helium 

Air 

Operating frequency (Hz) 

83 

0.2 

Pressure (MPa) 

2.59 

0.101 

Temperature of the hot end (K) 

923 

313 

Temperature of the cold end (K) 

353 

303 

t/ m ax of the hot end (m/s) 

3.7 

0.24 

Umax of the cold end (m/s) 

2.85 

0.24 

Kinematic viscosity of the hot end (nf/s) 

32.3xl0‘ 6 

15.9xl0' 6 

Kinematic viscosity of the cold end (m 2 /s) 

6.48xl0' 6 

15.9xl0' 6 

Reynolds number, Re ma x> of the hot end 

19.7 

74.5 

Reynolds number, Re ma x. of the cold end 

75.7 

74.5 

Valensi number, Va, of the hot end 

0.12 

0.47 

Valensi number, Va, of the cold end 

0.6 

0.47 


5.3.1 The LSMU Experimental Setup 

Figure 5.7 shows the set up of the experiments under unidirectional flow. At one side of the LSMU 
slices, the transition piece is connected with a fan by a 1 1.08 m (12 ft) long flexible tube and a 0.54 m 
(21 in.) long acrylic tube. At the other side, the transition piece connects the LSMU plates with a sector- 
of-an-annulus shaped opening. The transition piece is for transitioning from a round to a sector-of- 
annulus cross section. It consists of 9 layers with the sector-of-annulus shaped opening and one layer with 
the round opening, which are shown in figure 5.8. The thickness of one layer is 12.7 mm (0.5 in.). Screen 
material (not shown in fig. 5.8) is sandwiched between every two layers to help with the flow diffusion. 
The pressure drop across the LSMU layers is measured by a micro manometer. The hot-wire anemometer 
is used to measure the velocity of the outlet flow. The voltage readings of the anemometer are input to the 
multimeter and then collected by the computer. 
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Pressure tap 



Figure 5.7. — The schematic of the experimental setup of unidirectional flow test. 



Figure 5.8. — Two transition sections (without screens). 


5.3.2 Traverse of the Hot Wire Probe Over the Area of a Sector of an Annulus 

A program was written in the c programming language to traverse the hot-wire probe over the area of 
a sector of an annulus. The measurement grid is shown in figure 5.9. Dots show the locations at which 
velocity measurements were taken. In the radial direction, the increment is 9 mm . In the upper area, 
including areas la, lb, 2a and 2b, the angular increment is tt/ 60. In the lower area, the angular increment 
is jr/120. There are 348 grid points in total. The order in which the probe visited the various areas is la, 
2a, 2b, lb, 3 to 10. The probe visited the centroid, shown by a cross in figure 5.9, before and after each 
area was visited. This allowed a check for time variations. 
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Figure 5.9. — The measurement grid. 


5.3.3 Friction Factor 

LSMU layers'. The Darcy velocity in the LSMU layers (i.e., the approach velocity to the LSMU), was 
measured, with 8 plates or layers in the LSMU for this friction-factor measurement. The local velocity 
inside the channels, V, of the LSMU test section is calculated by dividing the Darcy velocity by the 
porosity. The Reynolds number is based on local velocity inside the channels and the hydraulic diameter 
of the channels, which is 4.87 mm. Also, the pressure drop, A p , across the LSMU test section was 
measured. The friction factor is obtained from: 
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(5.5) 


Apx(D h / L) 

■’ 1 

-pV 2 
2 

where L is the length, or sum of the LSMU layer thicknesses. The result will be compared with the 
friction factors of continuous and staggered channels of parallel plates, a random fiber matrix and a 
woven screen matrix in figures 5.1 1 and 5.12 (and later against CFD computational results; see sec. 6.0). 
The following paragraphs describe the origins of the friction-factor data to be compared with the 
measured LSMU friction factor. 

Continuous channels : First, the equivalent continuous channel geometry to the LSMU regenerator is 
described. The LSMU channel length varies from 54.5 to 32.7 mm for the 6 rib plate. The average length 
of 43.6 mm is chosen to calculate the aspect ratio for the continuous channel comparison case. Thus, the 
aspect ratio is 43.6/2.58 = 17. The hydraulic diameter of the equivalent continuous channels is 4.87 mm. 
For fully developed flow in a continuous channel, /*Re = 89.9 when the aspect ratio is 20;/ K Re = 96 for 
an infinite aspect ratio (Munson et al. 1994). The friction factor of fully developed flow of a 43.6 mm 
long by 2.58 mm wide continuous channel is calculated by interpolation as /= 89.9/Re. This is the 
continuous channel equation plotted in figures 5.1 1 and 5.12. 

Staggered fins or plates: The curve for the Fanning friction factor for a strip fin configuration from 
Kays and London (1964) is used to describe staggered fin channels. These “staggered plates” curves of 
figure 5.1 1 and 5.12 are taken from figure 5.10 which is copied from Kays and London (1964). The Darcy 
(or Darcy-Weisbach) friction factor is computed from the Fanning friction factor by multiplying by 4.0. 

Random fiber: From the Sage code manuals (Gedeon, 1999), the friction factor for a random fiber 
matrix is: 


/ = 192/Re+4.53Re -0067 (5.6) 

The local bulk mean velocity (not the Darcy velocity) is used to calculate the Reynolds number. 

Woven screens: From Sage (Gedeon, 1999), the friction factor for a woven screen matrix is: 

/ = 129/Re+2.91Re-°- 103 (5.7) 


Again, the local velocity (not the Darcy velocity) is used to calculate the Reynolds number. 

A comparison of friction- factor as a function of Reynolds number, or / versus Re, for the different 
geometries is shown in table 5.2. Figure 5.11 shows the Darcy-Weisbach friction factor versus Reynolds 
number for the different geometries. The Reynolds numbers are based upon the local, in-channel velocity 
and the hydraulic diameter is computed as . The curve for the continuous channel is based on 

4 wetted 

laminar flow. Figure 5.12 shows the same data but in a log plot. 
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Fin pitch = 19.82 per in. 

Plate spacing, b = 0.205 in. 

Splitter symmetrically located 
Fin length flow direction = 0. 125 in. 

Flow passage hydraulic diameter, 4 r h = 0.005049 ft 
Fin metal thickness = 0.004 in,, nickel 
Splitter metal thickness = 0.006 in. 

Total heat transfer area/volume between plates, 13 = 680 ft 2 /ft 3 
Fin area (including splitter )/total area = 0.841 

Figure 5.10. — (versus Re for staggered strip-fin. (Figure used with permission from 
Kays and London (1964).) 
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Staggered plates 
(Kays & London) 

» LSMU plates 

Woven screens 
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Re 

Figure 5.1 1 . — Darcy-Weisbach friction factor versus Reynolds number for 
different geometries. The Reynolds number is based on the local 
velocity and hydraulic diameter. 



— Continuous channel 

-*- Staggered plates 
(Kays & London) 

» LSMU plates 

Woven screens 
(Sage) 

Random fiber matrix 
(Sage) 


Figure 5.12. — Darcy-Weisbach friction factor versus Reynolds number for 
different geometries in log scale. The Reynolds number is based on the local 
velocity and hydraulic diameter. 
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TABLE 5.2— COMPARISON OF FRICTION-FACTOR, /VERSUS 
REYNOLDS NUMBER, Rc, FOR DIFFERENT GEOMETRIES 


Rc 

/ continuous 
channel 

/ staggered 
plates 

/LSMU 

layers 

/ random 
fiber 

/ woven screens 

40.00 

2.25 



8.34 

5.22 

50.00 

1.80 



7.33 

4.52 

60.00 

1.50 



6.64 

4.06 

70.00 

1.28 



6.15 

3.72 

131.55 

0.68 


1.29 

4.73 

2.74 

156.61 

0.57 


1.08 

4.45 

2.55 

150.00 

0.60 



4.52 

2.60 

199.13 

0.45 


0.86 

4.14 

2.33 

200.00 

0.45 

0.630 


4.14 

2.33 

268.30 

0.34 


0.69 

3.83 

2.12 

250.00 

0.36 

0.525 


3.90 

2.16 

340.68 

0.26 


0.56 

3.63 

1.97 

300.00 

0.30 

0.464 


3.73 

2.05 

408.97 

0.22 


0.53 

3.50 

1.88 

350.00 

0.26 

0.414 


3.61 

1.96 

477.13 

0.19 


0.48 

3.40 

1.81 

400.00 

0.22 

0.382 


3.51 

1.89 

450.00 

0.20 

0.349 


3.44 

1.84 

500.00 

0.18 

0.330 


3.37 

1.79 

600.00 

0.15 

0.294 


3.27 

1.72 

800.00 

0.11 

0.254 


3.13 

1.62 

1000.00 

0.09 

0.227 


3.04 

1.56 

1067.26 

0.08 


0.27 

3.02 

1.54 

1500.00 

0.06 

0.194 


2.90 

1.46 

1556.17 

0.06 


0.24 

2.89 

1.45 

2000.00 

0.04 

0.179 


2.82 

1.39 

2350.71 

0.04 


0.21 

2.77 

1.36 

3000.00 

0.03 

0.163 


2.71 

1.32 


5.3.4 Permeability and Inertial Coefficient 

The Darcy -Forchheimer equation (a steady-flow form of the 1-D momentum equation) can be written 
in terms of the empirical coefficients permeability, K, and inertial coefficient, C/ as 


A p 

L 


= il/ + 

K 


vz 


pU 2 


(5.8) 


Equation (5.8) could be rewritten as 


AutuU v 

LU K VZ 


(5.9) 


The Darcy velocity U can be calculated by 


C/ = ( p-V (5.10) 

where cp is the porosity. The porosity is 0.86 for the LSMU layers and is 0.9 for staggered plates (Kays 
and London, 1964). The porosity values for the woven screens and random fibers were both chosen to be 
0.9. The porosity dependence parts of equations (5.6) and (5.7) are given in Gedeon’s Sage manuals 
(Gedeon, 1999). The porosity of the continuous channels is chosen to be the same as that for the LSMU 
plates (0.86). 
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From the plot of ^ versus U, the intercept is taken to be . Permeability is calculated by 
dividing the dynamic viscosity by the intercept value. Figures 5.13 through 5.17 show the A 7(/ ^ versus 

U plots for continuous channels, staggered plates, LSMU layers, woven screens and random fibers, 
respectively. For continuous channels, there is no inertial effect (since the flow path contains no 

“obstacles”). The plot is a horizontal line and is taken to be 40.7. For staggered plates, is taken 

to be 434, found by extrapolation to zero velocity since no Darcy behavior was seen in the plot (no flat 

section in fig. 5.14). For the LSMU layers, is taken to be 77. Woven screens and random fiber 

matrices have a similar pattern to one another, an almost constant slope over the velocity range available. 
A fitting equation is generated for the lower-velocity region and is extrapolated to the U=0 m/s intercept 
for evaluation of K, as discussed above. The intercept is obtained from the fitting curve. For the woven 

screen and random fibers, is taken to be 36439 and 53201, respectively. Table 5.3 shows 

permeability for different geometries. 



♦ Continouous channels 


U (m/s) 

Figure 5.13 . — ApILIU versus l/for continuous channels. 



• Staggered plates 

Linear (Staggered 

plates) 


U (m/s) 

Figure 5.14 . — ApILIU versus U for staggered plates. 
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Figure 5.15 . — ApILIU versus L/for LSMU layers. 



U (m/s) 

Figure 5.16 . — ApILIU versus L/fora woven screen matrix. 
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U (m/s) 

Figure 5.17. — Ap/LIU versus U for random fibers. 


TABLE 5.3— COMPARISON OF PERMEABILITY 
FOR DIFFERENT GEOMETRIES 



Dh, 

mm 

K, 

m 2 

K/Dh 2 

Continuous channels 

4.872 

4.54E-07 

1/52.3 

Staggered plates* 

1.539 

4.26E-08 

1/55.6 

LSMU layers 

4.872 

2.40E-07 

1/98.98 

Woven screen matrix* 

0.2 

5.07E-10 

1/78.9 

Random fibers* 

0.2 

3.47E-10 

1/115.2 


* Extrapolated into Darcy region. 


The choice of the hydraulic diameter of the woven screen matrix and random fibers will not affect the 
value of K/Dh 2 , which is proven by the following. 

Equation (5.5) can be rewritten as 


Ap _ /• p-V 
V • L 2 D h 

Since Re = VDyf b and U = Ftp, 

Equation (5.11) can be written as 

Ap _ f ■ p • Re- u 
U L 2<p D h 2 


(5.11) 


(5.12) 


As shown by equation (5.9), when the flow velocity is small, 



term, 



jj is dominated by the viscous 


A p _ f • p • Re- u _ p 
U^L' 2 (p-D/t 2 
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(5.14) 


K _ 2<p 
D h 2 /'Re 

For random fiber matrix, from equation (5.6), 

/Re = 192 + 4.53 Re 0 - 933 (5.15) 

For a woven screen matrix, from equation (5.7), 

/Re = 129 + 2.91 Re 0 - 897 (5.16) 

Therefore, for small Reynolds numbers, the value of / Re is a constant, 192 and 129 for random 
fiber and woven screen matrix, respectively. Thus, KJDh 2 does not change with the selection of the 
hydraulic diameter of the matrix, provided that the porosity remains fixed. 

Equation (5.9) could be rewritten as 


A p 

LU 2 


KU 


y[K 


(5.17) 


The inertial coefficient, C f , can be evaluated accurately as the velocity goes toward infinity. However, 
due to the limitation of the equipment available, 6.5 m/s is the maximum mean velocity that was reached. 

Figures 5.18 and 5.19 show the ^ /2 versus U plot for LSMU plates and staggered plates, 


respectively. For LSMU plates, 


A/6.473 


-Jk 


is taken to be 34.6. For staggered plates, 




C 


/ 


H — — p is taken to be 77.5. The inertial coefficient values are 0.00939 for LSMU plates and 
K27.37 VZ 

0.01075 for staggered plates. 
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Figure 5.18. — Ap/LIU 2 versus L/for LSMU plates. 


NASA/CR— 2007-2 15006 


57 


350.000 n 



« Staggered plates 


U (m/s) 

Figure 5.19. — A pILIU 2 versus U for staggered plates. 


5.3.5 Friction Factor of Various LSMU-Plate Configurations 

Figure 5.21 shows the comparison of Darcy friction factor of various LSMU plate configurations, 
including 8 LSMU plates, 10 LSMU plates, 5 aligned LSMU plates, and 8 LSMU plates (double 
thickness). For the “aligned” test, five 6-rib LSMU plates are stacked together and tested under the 
unidirectional flow. Figure 5.20 shows a picture of 5 aligned LSMU plates. The fins are aligned 
throughout the entire area. One must take the view toward the bottom of figure 5.20 to see this. The total 
thickness of the 5 plates is 39.7 mm. The hydraulic diameter, d h , of the flow channel is 4.87 mm. The 
ratio of the length to the hydraulic diameter is 8.15. For laminar flow in a continuous channel, the ratio of 
the entrance length to the hydraulic diameter is 0.06*Re. Under current test conditions, the Reynolds 
number varied from 207 to 1618. Thus, the entrance length changes from 12.4(7/, to 97. \d h and the flow of 
the 5 aligned plates is in the developing regime. Figure 5.21 shows the aligned plates have lower friction 
factor values than those for the standard LSMU plate arrangement. This is because the flow through the 
aligned plates is continuous and has minimal flow separation whereas the LSMU plates under the 
standard arrangement woidd have wakes from trailing edges and separation on leading edges. 

A comparison can be made between the case where 10 LSMU plates are stacked under the standard 
configuration and the case where 8 LSMU plates are stacked similarly. The two cases compare very 
closely. The shorter assembly has only slightly larger friction factor values. This is an indication that the 
flow develops rapidly within the assembly, perhaps in the first 3 or 4 plates. The fitting equation for 
friction factor versus Reynolds number for the 8 LSMU plates is: 

/ = + o.i27 Re 001 (5.18) 

Re 

This equation fits the 1 0 LSMU plate data also. 


NASA/CR— 2007-2 15006 


58 



Figure 5.20. — A picture of 5 aligned LSMU plates. 



Figure 5.21. — Comparison of Darcy friction factors of various LSMU-plate configurations. 

To determine what the friction factor would be with the LSMU geometry but with plates that are 
twice as thick, we stacked two, 6-rib plates together and two 7-rib plates together, and then repeated, 
giving 4 groups with the eight LSMU plates (or the equivalent of 4 double -thickness plates). Figure 5.21 
shows that the friction factor is reduced from that with normal stacking with this new stacking order. The 
reason is that there are fewer flow redistributions from 6-rib geometry to 7-rib geometry, or the reverse, 
with this stacking order. 
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5.4 The Jet Penetration Study 

5.4.1 The Jet Penetration Study for the Round Jet Generator 

The geometry of the round jet generator is shown in figure 5.22. The fabricated round jet generator is 
shown in figure 5.23. Round holes are arranged in an equilateral triangle pattern. The hole diameter is 20 
mm and the center-to-center spacing is 40 mm. The jet generator is 30.5 cm (12 in.) long giving a hole 
L/D ratio of 15.2. 

5.4. 1.1 The Experimental Setup 

Since the diameter of the hot wire support tube is 4.57 mm, a spacer for hot wire insertion was 
fabricated that is thicker than 4.57 mm, smaller than, but comparable in size to, the thickness of an LSMU 
plate (7.95 mm or 5/16 in.). Insertion of the hot-wire probe is thus limited to the plenum between the jet 
generator and the matrix and velocity data are not taken between LSMU layers. Temperature data are 
used to document the thermal field in the matrix, as affected by the penetrating jet. For this measurement, 
a thermocouple wire, which is much thinner than the velocity probe, is passed through a much thinner 
spacer inserted between two adjacent LSMU layers in the test. The thermal effect of the jet flow within 
the regenerator is documented by temperature profiles from the thermocouple traversed in the radial 
direction (of the test section which is a sector of an annulus, see fig. 5.22) across the primary jet of the test 
section and its two neighboring jets on that radius (also shown in fig. 5.22). Traverses are taken at several 
axial locations to document the widening thermal effects of the jets as one moves away from the jet 
generator. The thermal field is due to the low temperature of the jet being dispersed within the matrix, 
hydrodynamically and thermally and, possibly, being turned radially away from its center by the matrix 
material. 


o o 



Figure 5.22. — The round jet generator and the 
plenum shape, a sector of an annulus. The 
center colored jet is the primary jet. 
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Figure 5.23. — The round jet generator. 
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Figure 5.24. — The schematic of the test facility. 


Figure 5.24 shows a schematic of the test facility. The main components of the facility are the 
oscillatory-flow generator, a cooler, two transition pieces (one on each end of the regenerator), a jet 
generator, ten LSMU plates, an electrical heating coil and an isolation duct. The cooler is a compact heat 
exchanger used to heat the passenger compartment of a car (called a heater core). One transition piece is 
put between the regenerator and the heating coil; the other transition piece is located between the jet 
generator and the cooler. The isolation duct is a long, open tube. It has active mixing to isolate the 
experiment from the room conditions. For the oscillatory flow generator, the stroke is 178 mm, the piston 
diameter is 216 mm and the frequency is 0.2 FIz, selected to match the Reynolds number and the Valensi 
number of the microfabricated regenerator in the pattern engine (79 and 0.53, respectively). There is a 
plenum between the jet generator and the LSMU plates which allows the hot-wire traverse for 
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documenting the periodicity of the jets along a radius through the centerlines of the holes highlighted in 
figure 5.22. The nominal thickness of the plenum, 8, is determined such that the axial-flow area, nd} 1 4, 
equates to the radial-flow area, nd c S . The diameter of the jet channel, d c , is 20 mm; thus, the nominal 

thickness for the spacer, 8, is 5 mm. This plenum allows a single hot-wire probe, with a support tube 
diameter of 4.57 mm, to be inserted into the plenum for taking velocity measurements. 

A spacer consisting of two 0.76 mm (0.030 in.) thick stainless-steel sheets (right and left of fig. 5.25) 
was inserted between two adjacent LSMU plates to allow the thermocouple wire used to take temperature 
profiles to pass through the test matrix. The opening of the spacer, which is the gap between two stainless 
steel sheets, is 0.51 mm (0.020 in.). 

Thermocouples of type E with a diameter of 76 pm (3 mils) are used for unsteady temperature 
measurements within the LSMU plate test section. The time constant of the thermocouple is 0.05 sec, 
which means the thermocouple can sufficiently quickly respond to changes in flow temperature. The 
thermocouple measuring the temperature within the LSMU matrix is mounted on a stepper motor driven 
rail so that it can be traversed inside the spacer and between two LSMU plates. The spacer can be moved 
to other axial locations within the LSMU matrix, allowing temperature documentation at various axial 
locations. The temperatures within the LSMU matrix, T(x, r, CA), are presented as functions of x, the 
distance in the axial direction; r, the distance along the centerline radius; and CA, the crank angle. One 
stationary thermocouple is located at one end of the jet generator and adjacent to the plenum. It is for 
measuring the cold end temperature, 7) (CA). Another stationary thermocouple is located at the end of the 



Slots for thermocouples 


Figure 5.25. — The spacer on the LSMU plates. 
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transition piece which is adjacent to the LSMU plates. It is for measuring the hot end temperature, 
7},(CA). At each location, these three temperatures are taken at a sampling frequency of 500 Hz for 50 
cycles. To eliminate temperature drift, a dimensionless temperature is calculated as: 


i()(x, r, CA) 


T(x,r,CA) - T c 

T h ~ T c 


(5.19) 


The cold end temperature, T c ( CA), is averaged over the portion of one cycle when the flow is passing 
from the jet generator to the LSMU regenerator plates. This gives an average temperature, T c , for each 
cycle. The hot end temperature, T h ( CA), is averaged over the portion of one cycle when the flow is 
passing from the heater to the LSMU regenerator plates. This gives an average temperature, T h , for one 
cycle. The dimensionless temperature <()(x, r, CA ) is calculated for each reading of the cycle. Averages of 
i))(x, r, CA ) are taken over an ensemble of 50 cycles. 

5.4.1.2 Jet-to-Jet Uniformity for the Round Jet Array 

To verify that the flow was uniformly distributed in the round jet generator under oscillatory flow 
conditions, the velocities within the plenum between the jet generator and the LSMU regenerator plates 
were measured. Results are given versus time within a oscillation cycle based upon an ensemble average 
of 50 cycles. Figure 5.22 shows the round jet generator and the plenum which is a sector of an annulus in 
shape. The hot-wire probe is driven by the stepper motor to move horizontally along a line which passes 
through the centerlines of the three holes highlighted in figure 5.22. Each velocity measurement is taken 
at a sampling frequency of 500 Hz for 50 cycles. Velocity profiles taken during the blowing half of the 
cycle, when the flow is passing from the jet generator to the LSMU plates, are shown in figure 5.26. The 
origin of the horizontal axis is the center of the center jet. Velocity profiles during the drawing half of the 
cycle, when the flow is passing from the LSMU plates back to the jet generator are shown in figure 5.27. 
Velocity profiles show that the jets from the three round channels shown are very similar to one another. 
This confirms that when the center jet is interrogated, the data are representative of data for flow from all 
interior jets. 



Location (mm) 


Figure 5.26. — Velocity profiles during the blowing half of the cycle, when the 
flow is passing from the jet generator to the LSMU plates. The origin of the 
horizontal axis is the center of the center jet. 
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Figure 5.27. — Velocity profiles during the drawing half of the cycle, when 
the flow is passing from the LSMU plates back to the jet generator. The 
origin of the horizontal axis is the center of the center jet. 


When the crank angle is 270°, the average velocity of the jet is 0.988 m/s. From the mass 
conservation equation for incompressible flow, the area-mean flow velocity within the round jet can be 
calculated: 


TC Dl 

U C A C = — ~U p (5.20) 

where A c is the open area of the jet generator, U c is the average velocity over the round jet generator, D p is 
piston diameter, and U p is the piston velocity. In the oscillatory flow generator, the piston moves in a 
sinusoidal fashion. The displacement, X p , can be calculated from the stroke and the frequency, 

X p = — - — cos(2 nft) (5.21) 

The piston velocity can be obtained by taking the first derivative of the piston displacement, 

Up — Xp — n /Stroke sin(27t / 1) (5.22) 

The opening area of the jet generator is 4308 mm 2 . The piston velocity is 

Up = 0.1 12sin(2rc ft) m/s 

The average velocity over the round jet generator is 

U c = 0.95sin(27i/f) m/s 

which matches very well with 0.988 m/s, which is the average velocity measured by the hot wire when 
the crank angle is 270°. 

5.4. 1.3 Some Important Parameters in the Round Jet Generator and the LSMU Plates 

The displacement of the fluid particle within the round jet also can be calculated from the mass 
conservation equation for incompressible flow, 
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(5.23) 


X C A C 


nDj 

4 


X 


P 


so, 


X c =757 cos {2n ft) mm 


The amplitude ratio, A R , is the fluid displacement during half a cycle divided by the tube length. For 
the round jet generator, which is 305 mm (12 in.) long, the amplitude ratio is 2.48. 

The maximum Reynolds number in the round jet generator is, 


Re r 


Uc max d 0.95x0.02 


= 1195 


v 15.9x 10 -6 

The Valensi number in the round jet generator is, 

d 2 xs 0.02 2 x2nx 0.2 

Va = = = 7.9 

4v 4 x 15.9 x 10 -6 

The flow velocity within the LSMU plates can be calculated from: 


U,.A r = 


nD: 


-u, 


(5.24) 


where A r is the open area of the LSMU plates. The flow velocity within the LSMU plates is: 

U r = 0.243 sin( 2nft ) m/s 

The displacement of the fluid particle within the LSMU plates can be calculated as, 


X y Ay ' 


7 xDi 


X r 


where 


X r = 194 cos(2n/t) mm 

For the LSMU of 10 plates, which is 79.4 mm long, the amplitude ratio is 2.44. 

5.4. 1.4 Jet Penetration of the Round Jet Generator 

Dimensionless temperatures at six axial locations have been measured. The six locations are: between 
the plenum on the jet generator side and the first LSMU plate, between plates 2 and 3, between plates 3 
and 4, between plates 5 and 6, between plates 8 and 9, and after plate 10. The dimensionless temperature 
profiles are shown in figures 5.28 through 5.33. 

During the blowing half of the cycle, when the flow is passing from the jet generator to the LSMU 
plates, the crank angle changes from 180° to 360°. Three cold jets, which are distinguished from the rest 
of the region, can be seen in figure 5.28. Figure 5.32 shows that the jet edges are nearly imperceptible. 

The jet penetration depth is about the thickness of 8 LSMU plates, which is 63.5 mm. The hydraulic 
diameter of the LSMU plates is 4.872 mm so the jet penetration depth is about 13 times the hydraulic 
diameter. A movie of the jet penetration was generated and can be obtained from the authors (contact 
tsimon@me.umn.edu.) In figures 5.32 and 5.33, one can see an effect of the boundary of the test section 
at the larger radius of figure 5.22 beginning to influence the temperature data for r > 35 mm. 
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Crank angle (degree) 



Location (mm) 

Figure 5.28. — Dimensionless temperature profiles between the plenum and the first 
LSMU plate. The origin of the horizontal axis is the center of the center jet. 
Neighboring jets are centered at -40 mm (the one at the larger radius of 
fig. 5.20) and 40 mm. 



Location (mm) 

Figure 5.29. — Dimensionless temperature profiles in the middle of the LSMU plates, 
between plate 2 and 3. The origin of the horizontal axis is the center of the center 
jet. 
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Location (mm) 


Figure 5.30. — Dimensionless temperature profiles in the middle of the LSMU plates, 
between plate 3 and 4. The origin of the horizontal axis is the center of the center jet. 



Location (mm) 


Figure 5.31. — Dimensionless temperature profiles in the middle of the LSMU plates, 
between plate 5 and 6. The origin of the horizontal axis is the center of the center 
jet. 
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Figure 5.32. — Dimensionless temperature profiles in the middle of the LSMU plates, 
between plate 8 and 9. The origin of the horizontal axis is the center of the center 
jet. 



Location (mm) 


Figure 5.33. — Dimensionless temperature profiles after 10 LSMU plates. The origin of 
the horizontal axis is the center of the center jet. 
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5.4.1. 5 


Jet Growth and the Fraction of Inactive Matrix Material 


The center jet’s edge, as the jet expands along the axial flow direction, is defined by assuming the 
edge occurs at that point at which the dimensionless temperature is the average of the maximum and the 
minimum dimensionless temperatures found in traversing across the jet— at a certain axial location and a 
certain crank angle. One way of representing this assumption about the jet width, b, is via the temperature 
expression: 


<K*> b / 2, CA) = j (<t> max (*, CA) + <j> min (x, CA)) (5.25) 

where ± b/2 represent the two edges of the jet along the radial direction, with the center of the jet at the 
origin. Throughout the blowing half cycle, this jet diameter remains almost invariant with crank angle. 

The dimensionless temperature at 270° cra nk angle is chosen to evaluate the jet growth. Figure 5.34 
shows the jet growth along the axial direction (Note that the jet edges are difficult to identify between 
plates 5 and 6 and the jet width there is extrapolated). 

Figure 5.35 shows the jet penetration in the matrix and the jet penetration depth x p . The fraction of 
inactive matrix material is the fraction of matrix material that is not participating fully in thermal 
exchange with the working medium over the jet penetration depth, x p . The fitting equation of figure 5.34 
can be used to get jet diameter, b{x), over ()<x/d h < 1 3 for the LSMU plates. For one jet, the corresponding 
matrix area is a hexagon with side length of 23.1 mm, which is shown in figure 5.36. The area of the 
hexagon is A r The fraction of inactive matrix material is calculated by 



A J~ 


n b(x) : 


dx 


A value of 47% is found for the LSMU plates with the round jet generator. 



Axial location (mm) 


Figure 5.34. — Jet width of the center jet at crank angle 270°. 


(5.26) 
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Figure 5.35. — Jet penetration. 



Figure 5.36. — Area assigned to each jet. 
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5.4.2 The Jet Penetration Study for the Slot Jet Generator 

Figure 5.37 shows the slot jet generator and the plenum which is a sector of an annulus in shape. The 
fabricated slot jet generator is shown in figure 5.38. The slot channels are separated by the fins in the slot 
jet generator. The channel width is 8.5 mm and the fin thickness is 23 mm. The jet generator is 30.5 cm 
(12 in.) long. The experimental setup is shown in figure 5.24. 

5.4.2. 1 Jet-to-Jet Uniformity of the Slot Jet Array 

To verify that the flow is uniformly distributed in the slot jet generator under oscillatory flow 
conditions, the velocities within the plenum between the jet generator and the LSMU regenerator plates 
are measured. Results are given versus time within an oscillation cycle based upon an ensemble average 
of 50 cycles. The hot-wire probe is driven by the stepper motor to move horizontally along a line which 
passes through the center of the jet generator, normal to the slots (as indicated in fig. 5.37). 

The velocity measurement is taken at a sampling frequency of 500 Hz for 50 cycles. Velocity profiles 
taken during the blowing half of the cycle, when the flow is passing from the jet generator to the LSMU 
plates, are shown in figure 5.39. The origin of the horizontal axis is the center of the center jet. Velocity 
profiles show that the jets from the three slot channels shown are very similar to one another. Velocity 
profiles during the drawing half of the cycle, when the flow is passing from the LSMU plates back to the 
jet generator, are shown in figure 5.40. These velocity profiles show that the velocities in the center are 
slightly lower than the velocities which are distant from the center. Figure 5.41 shows the geometry of the 

6- rib LSMU plate and the traversing route of the hot-wire probe. The flow in the center is very close to 
the rib of the 6-rib LSMU plate, which decreases the flow velocity. Figure 5.42 shows the geometry of the 

7- rib LSMU plate and the traversing route of the hot wire probe. 



Figure 5.37. — The slot jet generator and the plenum shape, a sector of an annulus. 
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Figure 5.38. — The slot jet generator. 
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Figure 5.39. — Velocity profiles during the blowing half of the cycle, when the flow 
is passing from the jet generator to the LSMU plates. The origin of the 
horizontal axis is the center of the center jet. 
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Figure 5.40. — Velocity profiles during the drawing half of the cycle, when the flow is passing 
from the LSMU plates back to the jet generator. The origin of the horizontal axis is the 
center of the center jet. 



Figure 5.41 . — Geometry of the 6-rib LSMU plate. 
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Figure 5.42. — Geometry of the 7-rib LSMU plate. 

When the crank angle is 270°, the average velocity of the jet is 0.783 m/s. From the mass 
conservation equation for incompressible flow, the flow velocity within the slot jet can be calculated: 

7 iDl 

UhA h - ~^~U p (5.27) 

where A h is the open area of the slot jet generator, U h is the average velocity over the slot jet generator, 

D p , is piston diameter and U p is the piston velocity. 

In the oscillatory-flow generator, the piston moves in a sinusoidal fashion. The displacement, X p , can 
be calculated from the stroke and the frequency, 


.. Stroke .. . . 
X p = cos(27t Jt) 


(5.28) 


The piston velocity can be obtained by taking the first derivative of the piston displacement, 


Up - X p =7 if Stroke [sin (2n f t )] (5.29) 

The open area of the slot jet generator is 5278 mm 2 . The piston velocity is 

Up = 0.1 12sin(27r/t) m/s 
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The average velocity over the slot jet generator is 

U c = 0.776 sin (inft) m/s 

which matches very well with 0.783 m/s, which is the average velocity measured by hot wire when the 
crank angle is 270°. 

5.4.2.2 Some Important Parameters in the Slot Jet Generator 

The displacement of the fluid particle within the slot jet also can be calculated from the mass 
conservation equation for incompressible flow, 


X h A h 


7 iD: 


■I, 


(5.30) 


so, 


X /, =618 cos(2tc / 1 ) mm . 


The amplitude ratio, A R , is the fluid displacement during half a cycle divided by the tube length. For 
the slot jet generator, which is 304.8 mm (12 in.) long, the amplitude ratio is 2.03. 

The maximum Reynolds number in the slot jet generator is, 


Re 


max 


u h . 


max 


d 


V 


0.776x0.017 

15.9xl0- 6 


= 830 


The Valensi number in the slot jet generator is, 

,, d 2 co 0.017 2 x 2rtx 0.2 r „ 

Va = = — = 5.7 

4v 4xl5.9xl0- 6 

5.4.2.3 Jet Penetration of the Slot Jet Generator 

Dimensionless temperatures at five axial locations were measured. The five locations were: between 
the plenum on the jet generator side and the first LSMU plate, between plates 3 and 4, between plates 5 
and 6, between plates 6 and 7, and between plates 8 and 9. The dimensionless temperature profiles are 
shown in figures 5.43 through 5.47. A movie of the jet penetration was generated and can be obtained 
from the authors (contact tsimon@me.umn.edu.) 

During the blowing half of the cycle, when the flow is passing from the jet generator to the LSMU 
plates, the crank angle changes from 180° to 360°. Three cold jets, which are distinguished from the rest 
of the region, can be seen in figure 5.43 (between the plenum and the first LSMU plate). Figure 5.46 
shows that the jet edges are nearly imperceptible between plates 6 and 7. The jet penetration depth is 
about the thickness of 6 LSMU plates, which is 47.6 mm. The hydraulic diameter of the LSMU plates is 
4.872 mm so the jet penetration depth is about 10 times the hydraulic diameter. In figures 5.44 and 5.45, 
one can see the effect of the rib of the 6-rib LSMU plate, which causes the center jet to have slightly 
higher temperatures during the blowing half cycle. 
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Figure 5.43. — Dimensionless temperature profiles between the plenum and the 
first LSMU plate. The origin of the horizontal axis is the center of the center jet. 
Neighboring jets are centered at -31 .5 mm (left side of fig. 5.37) and 31 .5 mm. 
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Figure 5.44. — Dimensionless temperature profiles in the middle of the LSMU plates, between plates 
3 and 4. The origin of the horizontal axis is the center of the center jet. 
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Figure 5.45. — Dimensionless temperature profiles in the middle of the LSMU plates, 


between plates 5 and 6. The origin of the horizontal axis is the center of the center jet. 
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Figure 5.46. — Dimensionless temperature profiles in the middle of the LSMU plates, 
between plates 6 and 7. The origin of the horizontal axis is the center of the center jet. 
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Figure 5.47. — Dimensionless temperature profiles in the middle of the LSMU plates, between plates 
8 and 9. The origin of the horizontal axis is the center of the center jet. 


5.4.2.4 Jet Growth and the Fraction of Inactive Matrix Material 

Throughout the blowing half cycle, the slot jet width remains almost invariant with cra nk angle. The 
dimensionless temperature at 270° crank angle is chosen to evaluate the jet growth. Figure 5.48 shows the 
jet growth. The jet edges are difficult to identify between plates 6 and 7 and the jet width there is found by 
extrapolation. 

The fraction of inactive matrix material is calculated by 

F-Tp-Wlkfc (5.31) 

x„ J SL 
p o 

where b{x) is the jet width, x p is the jet penetration depth, S is the jet center-to-center spacing (31.5 mm), 
and L is the jet length. A value of 69% is found for the LSMU plates with the slot jet generator. This 
compares to 47% for the round jets entering the LSMU regenerator and 55% for the round jets entering 
the 90% porous screen regenerator studied by Niu et al. 2003. 

5.5 Unsteady Heat Transfer Measurements 

Unsteady heat transfer with the LSMU plates was investigated. The LSMU dynamically simulates the 
microfabricated regenerator plates of the segmented-involute-foil regenerator. 
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5.5.1 Embedded Thermocouple 

Figure 5.49 shows a sketch of one embedded thermocouple. The drilled hole is 0.30 mm (0.012 in.) 
diameter and 2 mm (0.080 in.) deep. Figure 5.50 shows a picture of one of the embedded thermocouple 
installations. 



Slot jet 
■ Round jet 


Axial location (mm) 

Figure 5.48. — Jet width of slot jet and round jet at crank angle 270°. 




Figure 5.50. — Picture of one embedded thermocouple. 
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5.5.2 Experimental Procedure 

There are three thermocouples mounted in the 6-rib plate and three thermocouples mounted in the 7- 
rib plate. The three thermocouple locations are labeled 1 to 3 along the radial direction pointing to the 
center of the arc, as shown in figure 5.51. For the case presented in this report, ten plates are stacked in 
the design order. The 6-rib plate with thermocouples is plate 5 and the 7-rib plate with thermocouples is 
plate 6. One thermocouple is traversed between plates 5 and 6 to measure the air temperature, as shown in 
figure 5.52. The thermocouple junctions in plate 5 are near the traversing thermocouple, while the 
thermocouple junctions in plate 6 are more distant from the traversing thermocouple. For every embedded 
thermocouple location, air temperatures are measured at 5 locations on the same radial line as that on 
which the embedded thermocouples reside: -2 mm, -1 mm, 0, 1 mm, and 2 mm away from the embedded 
thermocouple radial location, labeled “1” through “5,” respectively. Two runs with different warm-up 
times, both sufficiently large, were conducted. Table 5.4 shows the case names of the two runs. Case B13 
is chosen to show the data processing. 



Figure 5.51 . — Locations of the embedded thermocouples. 


Figure 5.52. — Temperature measurement 
locations for the case presented. 
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TABLE 5.4.— CASE NAMES OF THE TWO RUNS 



Case 

Date 

Number of plates 

Tra. t.c. between plate 

Embedded t.c. loc. 

Tra. t.c. loc. 

1 

All 

8/11 

10 

5 and 6 

1 

1 

2 

A12 

8/11 

10 

5 and 6 

1 

2 

3 

A13 

8/11 

10 

5 and 6 

1 

3 

4 

A14 

8/11 

10 

5 and 6 

1 

4 

5 

A15 

8/11 

10 

5 and 6 

1 

5 

6 

A21 

8/11 

10 

5 and 6 

2 

1 

7 

All 

8/11 

10 

5 and 6 

2 

2 

8 

All 

8/11 

10 

5 and 6 

2 

3 

9 

A24 

8/11 

10 

5 and 6 

2 

4 

10 

A25 

8/11 

10 

5 and 6 

2 

5 

11 

A31 

8/11 

10 

5 and 6 

3 

1 

12 

A32 

8/11 

10 

5 and 6 

3 

2 

13 

A3 3 

8/11 

10 

5 and 6 

3 

3 

14 

A34 

8/11 

10 

5 and 6 

3 

4 

15 

A35 

8/11 

10 

5 and 6 

3 

5 

16 

B 1 1 

8/16 

10 

5 and 6 

1 

1 

17 

B12 

8/16 

10 

5 and 6 

1 

2 

18 

B13 

8/16 

10 

5 and 6 

1 

3 

19 

B14 

8/16 

10 

5 and 6 

1 

4 

20 

B15 

8/16 

10 

5 and 6 

1 

5 

21 

B21 

8/16 

10 

5 and 6 

2 

1 

22 

B22 

8/16 

10 

5 and 6 

2 

2 

23 

B23 

8/16 

10 

5 and 6 

2 

3 

24 

B24 

8/16 

10 

5 and 6 

2 

4 

25 

B25 

8/16 

10 

5 and 6 

2 

5 

26 

B31 

8/16 

10 

5 and 6 

3 

1 

27 

B32 

8/16 

10 

5 and 6 

3 

2 

28 

B33 

8/16 

10 

5 and 6 

3 

3 

29 

B34 

8/16 

10 

5 and 6 

3 

4 

30 

B35 

8/16 

10 

5 and 6 

3 

5 


After at least 5 hr of warm up time, data collection begins: 

Step 1 : Hot end and cold end plenum temperatures of the LSMU test setup are collected for 20 cycles. 

Step 2: Solid temperatures at location 1 in the 6-rib plate and the 7-rib plate, and air temperature around 
location 1 are collected simultaneously. Data are measured over 50 cycles. 

Step 3: Hot end and cold end plenum temperatures of the LSMU test setup are collected for 20 cycles. 

Step 4: Solid temperatures at location 2 in the 6-rib plate and the7-rib plate, and air temperature around 
location 2 are collected simultaneously. Data are measured over 50 cycles. 

Step 5: Hot end and cold end plenum temperatures of the LSMU test setup are collected for 20 cycles. 

Step 6: Solid temperatures at location 3 in the 6-rib plate and the 7-rib plate, and air temperature around 
location 3 are collected simultaneously. Data are measured over 50 cycles. 

Step 7: Hot end and cold end plenum temperatures of the LSMU test setup are collected for 20 cycles. 
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5.5.3 LSMU Unsteady-Heat-Transfer-Measurement Results 

Assuming the axial temperature distribution of the fin is linear, one can perform an energy balance of 
the plate fin in the vicinity of the embedded thermocouple: 


h(x,r,t)A s (Tf(x,r,t)-T s (x,r,t )) = mC (5.31) 

ot 

where h is the convective heat transfer coefficient, A s is the surface area of the plate, 7) is the air 
temperature, T s is the temperature of the plate, m is the mass of the plate, and C is the specific heat of the 
plate material. Equation (5.31) becomes: 


h(x, r, t)(T f (x, r, t ) - T s (x, r,t)) = pC?~ (5.32) 

2 ot 

where 5 is the plate thickness and p is the density of plate material. The convective heat transfer 
coefficient can be calculated from the measured temperature differences between the air and the plate and 
temporal gradients of the plate temperature. The Nusselt number can be obtained from: 

Nu = h^JL (5.33) 

k 

where d h is the hydraulic diameter of the channel (4.87 mm) and k is the thermal conductivity of the air. 

Figure 5.53 shows the air temperature at location 3 and the solid temperatures of the 6-rib plate and 
the 7-rib plate at location 1 in case B13. Consider the axial difference between the traversing 
thermocouple junction and the embedded thermocouple junction in the nearest plate (the 6-rib plate). 
Because there is an axial gradient in the system, the air temperature at the location of the air thermocouple 
is not the air temperature at the axial location of the plate thermocouple junction. A small correction is 
made. It is assumed that the axial temperature gradient of the air can be obtained from the temperature 
difference between the 6-rib and 7-rib plates. This temperature gradient is used to shift the air temperature 
from the air thermocouple location 1.3 mm (0.05 in.) to the embedded thermocouple location (6-rib 
plate). Figure 5.54 shows the temperature difference between the air and solid after the shift. This 
temperature difference profile is not balanced, which means the value of the peak does not match the 
valley. Figure 5.55 shows the comparison of the Nusselt number of current experiment with the 
correlation from the NASA/Sunpower oscillating-flow test rig. Recall that similar measurements by Niu 
et al. (2003), but in a wire-screen matrix, showed a similar plot of Nusselt number versus cycle position. 
The following features were noted: The heat flux computed from the solid temporal gradient is zero when 
the temperature difference is not, creating a zero Nusselt number. The temperature difference becomes 
zero when the heat flux is not, creating an infinite Nusselt number. When Niu et al. (2003) compared the 
measured results to correlation results computed by assuming quasi-steady behavior, the comparison was 
close only when the fluid velocity was near the peak value or during the deceleration part of the cycle. We 
expected similar behavior here. 

Figure 5.56 shows the comparison of the heat flux from this LSMU experiment with the correlation 
from the NASA/Sunpower oscillating-flow test rig. The lack of symmetry of Nusselt number and heat 
flux of the current experiment results from the lack of symmetry of the temperature difference profile. A 
balanced temperature difference profile can be generated by a small shift, moving the curve of figure 5.54 
vertically 0.04 °C. This is within the uncertainty in measuring a temperature within this system. This 
gives the symmetric curve we expect (since the measurement is in the axial center of the LSMU plates). 
Figure 5.57 shows the temperature difference between the air and solid based after this shift. Figures 5.58 
and 5.59 show the Nusselt number comparison and heat flux comparison after this shift is made. 
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Figure 5.53. — Air temperature at location 3 (V air ), solid temperatures of the 
6-rib plate (7 s6 ), and the 7-rib plate (T s 7), at location 1 in case B13. 



Figure 5.54. — The temperature difference between the air and solid. 
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Figure 5.55. — Comparison of the Nusselt number of the current LSMU 

experiment with the correlation from the NASA/Sunpower oscillating-flow 
test rig. 



Crank angle (degree) 


Figure 5.56. — Comparison of the heat flux of the current LSMU experiment 
with that given by the correlation from NASA/Sunpower test rig and the 
measured temperature difference. 
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Figure 5.57. — The temperature difference between the air and solid, with the 
shifted air temperature. 



Figure 5.58. — Comparison of the Nusselt number of current LSMU 

experiment (with the shifted air temperature) with the correlation from 
NASA/Sunpower oscillating-flow test rig, equation (3.4). 
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Figure 5.59. — Comparison of the heat flux of the current LSMU experiment 
with the correlation from the NASA/Sunpower oscillating-flow test rig and 
the measured temperature difference with the shifted air temperature. 


In the following, only the original temperature shift is applied (to effectively move the air 
thermocouple to the axial location of the metal thermocouple). The Nusselt numbers of the 6 cases for 
which the traversing thermocouple is at location 3 (closest to the embedded thermocouple) are calculated 
for cycle positions 120° and 300°. The results are shown in table 5.5. The average value is 7.14 and the 
RMS is 1.07. 


TABLE 5.5— NUSSELT NUMBER FOR DIFFERENT CASES 


Case 

Crank angle 120° 

Crank angle 300° 

A13 

7.7302 

7.5696 

A23 

8.2965 

6.545 

A3 3 

6.5366 

5.3525 

B13 

7.3799 

8.6731 

B23 

7.2222 

7.9134 

B33 

5.1717 

7.2729 


For the following, the second temperature shift is applied (to make the temperature difference plot 
symmetric, essentially to get the mean air temperature and the mean metal temperature equal to one 
another). The Nusselt numbers of the 6 cases for which the traversing thermocouple is at location 3 
(closest to the embedded thermocouple), are calculated at cycle positions 120° and 300°. The results are 
shown in table 5.6. The average value is 7.07 and the RMS is 0.86. We note that at a cycle position of 
300° the local velocity in our test is 0.21 m/sec. With this, the correlation of Gedeon (eq. (3.4)) gives a 
Nusselt number of 9.1. The average value in table 5.6 is 7.07, which is about 22% lower than the Gedeon 
value. The difference is 2.4 standard deviations so we expect it to be significant. We note that the 
microfabricated regenerator had roughness at the entrance of each channel due to EDM debris whereas 
our LSMU did not. Roughness would tend to enhance heat transfer. 
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TABLE 5.6.— NUSSELT NUMBER FOR DIFFERENT 


CASES BASED ON THE NEW SHIFT 


Case 

Crank angle 120° 

Crank angle 300° 

A13 

7.7302 

7.5696 

A23 

7.3766 

7.2917 

A33 

6.0522 

5.7459 

B13 

8.2363 

7.6841 

B23 

7.7836 

7.3048 

B33 

6.1025 

5.9403 


6.0 Analysis Tools and CFD Results 
(Cleveland State University (CSU)) 

6.1 Computational Fluid Dynamics (CFD) Summary 

The microfabricated, segmented-involute-foil regenerator was numerically investigated utilizing 
commercial CFD software under both steady-state and oscillatory-flow conditions. The geometry consists 
of a stack of disks (segments) with each disk containing involute-shaped, micron-range channels, with 
channel flow direction perpendicular to the plane of the disk. The lateral orientation of the channels 
alternates from disk to disk in the flow direction. Two-dimensional (2-D) and three-dimensional (3-D) 
simulations were carried out. Steady-state simulations were performed for Reynolds numbers from 50 up 
to 2000 based on the channel hydraulic diameter and the mean flow-field velocity. 

The results of this CFD research have been validated by comparing the CFD data with the literature 
and experimental correlations obtained at UMN (LSMU, large-scale geometry) and at Sunpower (actual- 
scale geometry). For the oscillatory-flow simulations in 2-D and 3-D a base case was chosen using helium 
as a working fluid, stainless steel for the solid material, 3 10 K at the hot end, 293 K at the cold end, a 
maximum Reynolds number, Re max , of 50 and a Valensi number, Re^of 0.229. The effects of changing: 

1) the oscillation amplitude and frequency, 2) the thermal contact resistance between layers, and 3) the 
solid material, on the total-regenerator heat loss (convection and conduction) were documented. These 
results are expected to be useful for further development of involute-foil regenerators for Stirling engines. 

6.2 CFD Introduction 

The Fluent CFD commercial code was used for 2-D and 3-D, steady and unsteady, fluid flow and heat 
transfer simulations of a microfabricated involute-foil regenerator for a Stirling engine. The knowledge 
gained enabled fundamental understanding of how fluid flow and heat transfer takes place inside the 
segmented-involute-foil flow paths. It also helped provide support for physical testing (large-scale and 
actual-size); comparison of the microscale CFD results with the test results and Sage 1-D code 
simulations, provided additional insight for making decisions about the involute-foil design details. 

In a Stirling engine the regenerator is one of the most important parts. In search of an improved 
regenerator, a novel segmented-involute-foil design was proposed. It consists of two types of alternately 
stacked disks with microfabricated channels. Figure 6.1 (similar to fig. 3.1 shown earlier) shows a 
progressively exploded flow direction view of one of the involute-foil regenerator disks. On the second 
zoom the channels can be seen more clearly. The disk portrayed in figure 6. 1 has six ribs. There is a 
second type of disk that has seven ribs. This enables the staggering of the ribs in the stack, to reduce axial 
conduction and improve axial channeling of the flow. 

Figure 6.2 shows a two-layer view down the flow path while figure 6.3 shows a 3-D view in which 
one can distinguish the alternate orientations of four layers of the involute-foil flow channels. Both 
figures show the alternating-layer ribs and indicate how each layer of ribs was offset from its two adjacent 
layers. Figure 6.4 (similar to fig. 3.2 shown earlier) shows a 3-D view of only one channel of one layer 
with dimension labels. It can be seen in this figure that the walls of the channel are curved with an 
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involute-foil profile. Table 6.1 shows involute-foil channel dimensions and the segment or layer length 
(thickness, in the table). Note that the channel “width,” W, given in table 6. 1 is a nominal value; as can be 
seen in the middle view of figure 6.1 , W varies somewhat from the inner to the outer ring of each layer of 
the actual-size involute-foil segments (this variation of U W” is also true for the large-scale segmented- 
involute-foils of the LSMU — as can be seen in figs. 5.41 and 5.42). 



Figure 6.1. — Exploded view of the microfabricated disk. 



rant layer rib 


Back layer rib 
~ midway 
between two 
front layer ribs 


Front layer rib 


Channel wall 


Figure 6.2. — Frontal view of two layers of microfabricated segmented-involute-foil channels. 
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Figure 6.3. — 3-D view of four layers of microfabricated segmented-involute-foil channels. 



Figure 6.4. — Segmented-involute-foil channel, of one layer. 


TABLE 6.1— INVOLUTE-FOIL CHANNEL DIMENSIONS 


Dimension 

Unit 

Value 

Gap, g 

pm, 10' 6 m 

86 

Gap+wall, s 

Jim 

100 

Wall thickness, s-g 

\im 

14 

Channel width, W 

Jim 

1000 

Disk (layer, segment) thickness, Lc 

jim 

265 

Porosity 


0.838 

Hydraulic diameter, Dh, 4A/P 

Jim 

162 


This design has several potential advantages compared to existing designs (such as random-fiber and 
wire-screen matrices) as demonstrated by Sage and Finite Element Analysis (FEA) analyses done in the 
Phase 1 part of this contractual effort — and the actual microfabricated hardware and test results from 
Phase II : 
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1) Improved ratio of heat transfer to pressure drop, i.e., high Figure of Merit, 

2) Low axial conduction due to minimal contact area between disks, 

3) Better reproducibility and control over geometric parameters, 

4) High structural integrity and durability. 

It was a challenge to analyze this geometry for fluid flow and heat transfer via CFD so as to facilitate 
the comparison between this geometry and other regenerator geometries. The first difficulty is the 
complex 3-D geometry, as depicted in figure 6.2. The second is the oscillatory nature of the flow as it 
occurs in the Stirling engine. 

In this section a description will be given for the modeling set-up used to analyze the problem. This 
includes the grid generation for the computational domain and the matrix under investigation. Results of 
the numerical investigation for the different parameters of the proposed problem will be presented 
followed by the conclusions. 

6.3 Establishing the Computational Domain 

It was decided early on that it would not be feasible from a microscopic computational point of 
view to model the whole regenerator. Therefore it was necessary to look for simplifications, which 
usually come in the form of symmetries and boundary-condition approximations. 

6.3.1 Radial Direction Periodicity 

One simplification comes from recognizing the periodicity in the radial direction that comes from the 
concentric arrangement of several rings of channels. The flow through the whole disk can be 
approximated by the flow through just one ring of channels situated halfway between the OD and the ID 
of the annulus (of the disk). A seven- fold reduction of the computational domain can be achieved with 
large savings of computational resources. In fact, without this reduction of the domain, the modeling of 
the geometry is not feasible. Furthermore, this simplification enables the next simplification. 

6.3.2 Angular Direction Periodicity 

This simplification comes from recognizing that there exists a sector of the ring of channels that, 
when repeated in the angular direction, reassembles the full ring. This is sometimes called circular 
symmetry. For the case studied the sector is about 8.87° and it reduces the domain by about 40 times 
(about 40 sectors make a full ring). The following figures show how these two periodicities are employed. 
Figure 6.5 shows a ring of channels selected from the middle of the microfabricated disk. 

By inspecting the ring of channels one can select for modeling only one sector and employ periodic 
boundary conditions as shown in figure 6.6. 

Figure 6.7 shows an enlarged area from the middle of the sector of figure 6.6. In figure 6.7, one can 
see the channel walls, the angle formed between channel walls (8 1 ° in this middle ring) in two successive 
layers and how the involute-foil profile of the wall deviates from a flat wall (by about 2°). This figure 
suggests possible further simplification of the geometry by approximating the involute-foil profile with a 
straight line and the angle between successive layer walls with a right angle (this will be discussed later). 
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Periodic sector (8.87°) 
(circular symmetry) 


Enlarged in Figure 6.7 


Figure 6.5. — Ring of channels in the “radial middle” of the 
microfabricated disk-annulus, and a sector of that ring. 


Figure 6.6. — Periodic sector from figure 6.5, showing two layers, with computational grids 
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Figure 6.7. — Enlarged area in the middle of the periodic sector. 


6.3.3 Flow Direction Periodicity 

This simplification comes from recognizing that the regenerator is a stack of just two types of 
alternating layers. Thus, the repeating unit is comprised of two layers. One can use the flow output of one 
repeating unit as the input to the next one, and so on. Figure 6.8 shows an isometric view of three 
successive layers (or more precisely, four boundaries of these three layers). Because the frontal area 
occupied by the circular ribs is very small in comparison to the rest of the frontal area, a simplification 
was made to line up the ribs from disk to disk. In the real geometry (see figs. 6.2 and 6.3), the ribs are 
staggered from disk to disk. The orientation (crossing angle) of the channel walls was not significantly 
altered by the alignment of the ribs. As shown in figure 6.8, the ribs of two successive layers are aligned 
but the channel walls are still approximately perpendicular from layer to layer. Aligning the ribs enables a 
bounded domain in the radial direction for both layers that forms a repeating unit. 

As mentioned above, the minimum thickness of the repeating unit has to be the thickness of two 
layers. However, the interface between two layers is a geometric discontinuity. The exit (velocities and 
temperatures profiles) of one repeating unit would be used as a boundary inlet to next one. It is better to 
have no geometric discontinuities at boundary inlets and outlets. Therefore the selected repeating unit 
consisted of half the thickness of one layer, followed by a full- thickness layer, and ending in another half 
thickness of the next layer. So a half-layer thickness was used at the entry and exit. Figure 6.9 shows this 
arrangement. 

6.3.4 Computational Domains for Oscillatory-Flow Simulations 

Flow direction periodicity works only for steady-state modeling. The transient simulation for the case 
studied requires oscillatory (alternating- direction, zero-mean) flow, requiring a stack of several layers to 
be included in the domain. A minimum of six layers was determined to be adequate to capture the 
oscillatory-flow phenomena. However even if radial and angular periodicities are employed, the grid size 
would still be too large for the available computation capability. Further simplifications, as indicated in 
figure 6. 10, must be used. If the foil-crossing angle from one layer to the next is approximated to 90° 
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( instead of 81), the involute-foil profile is approximated as straight and the round ends of the channels are 
neglected, then one can build a manageable grid. This is expected to capture most of the 3-D oscillatory- 
flow phenomena of the microfabricated design. Figure 6.10 shows such a computational domain. 


Inlet 




Figure 6.9. — 3-D involute-foil-layers computational domain, entry unit and repeating unit. 
The half-full-half computational domain unit can be repeated periodically until the 
required stack height is achieved. 
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Figure 6.10. — 3-D straight-channel-layers computational domain, for 6 layers. 


Table 6.2 shows dimensions of the computational domain shown in figure 6. 10. In order to maintain 
the same hydraulic diameter as the actual involute-foil geometry, the gap was adjusted to 81 pm, from 84 
(see section 3.6.2. 1). In order to maintain the same spacing, the wall thickness was adjusted to 19 pm, 
from 16. Another adjustment was made to, the layer thickness which was decreased from 265 to 250 pm. 
This was done to better match the actual disks that were fabricated for experimental testing, since the 
original design called for 265 pm thick disks, but the manufacturer (Mezzo) fabricated 250 pm thick 
disks. The resulting porosity, as shown in table 6.2, was 0.81 instead of the actual regenerator’s 0.84. 


TABLE 6.2— 3-D STRAIGHT CHANNEL 
COMPUTATIONAL DOMAIN DIMENSIONS 


Dimension 

Unit 

Value 

Gap, g 

pm, 10 6 m 

81 

Gap+wall, s 

jam 

100 

Wall thickness, s-g 

jam 

19 

Channel width, W 

jam 

300, symmetry 

Disk (layer) thickness, Lc 

jam 

250 

Porosity 


0.81 

hydraulic diameter, Dh, 2g 

jam 

162 


6.3.5 Two-Dimensional Computational Domains 

Further simplifications were also made, for part of the CFD study, by using a 2-D computational 
domain. For the cases studied the 2-D domain consisted of a single parallel-plate channel with 6 
successive sections. There was no variation in flow geometry upon exiting one section and entering the 
next. Flowever, by changing the solid-interface settings one could set various values for the thermal 
contact resistance (TCR) between sections. This was expected to capture the interruption in the wall 
thermal conduction that is obtained by alternating the orientation of the channel walls (from one layer to 
the next) in the 3-D domain. Figure 6.11 shows such a 2-D computational domain. This 2-D domain 
allowed for quick parametric studies and finding trends that could be later confirmed in 3-D with fewer 
runs. 
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Thermal contact resistance between layers 



As presented above, three computational domains, two 3-D and one 2-D, were identified for dealing 
with the geometry to be studied. They represent different levels of compromise between the actual 
problem and the resources required to model the problem. They are at the first level of a study matrix. 
The next levels are formed by considering different boundary conditions to be imposed on these 
computational domains. 

6.4 Material Properties and Boundary Conditions 
6.4.1 Material Properties 

For the fluid, helium gas was used at an operating pressure of 2,500,000 Pa. Table 6.3 shows the 
properties used for helium. 

As for the solid, two materials were used: stainless steel and pure nickel with properties as shown in 
tables 6.4 and 6.5, respectively. 


TABLE 6.3— HELIUM PROPERTIES USED IN THIS STUDY 


Property 

Units 

Method | Values j | j 

Density 

kg/m 3 

incompressible ideal-gas 

Cp 

J/kg-K 

constant 

| 5193 | 




T u coef. 

T 1 coef. 

T 2 coef. 

T 3 coef. 

Thermal conductivity 

W/m-K 

polynomial 

0.01998 

0.000509 

-2.61E-07 

8.53E-1 1 




Cl 

C2 



Viscosity 

kg/m-s 

Sutherland Law 

1.46E-06 

79.96 



Molecular weight 

kg/kg-mol 

constant 

| 4.0026 | 


TABLE 6.4— STAINLESS-STEEL PROPERTIES 


Property 

Units 

Method 

Values 




Density 

kg/m 3 

constant 

8030 







T° coef. 

T 1 coef. 

T 3 coef. 

T 3 coef. 

Cp 

J/kg-K 

polynomial 

148.86 

1.5139 

-0.0018 

7.35E-07 

Thermal conductivity 

W/m-K 

polynomial 

6.182 

0.03505 

-2.63E-05 

1.02E-08 


TABLE 6.5— NICKEL PROPERTIES 


Property 

Units 

Method 

Values 

Density 

kg/m 3 

constant 

8900 

Cp 

J/kg-k 

constant 

460.6 

Thennal conductivity 

W/m-k 

constant 

91.74 
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6.4.2 Dimensionless Quantities, Correlations and Boundary Conditions 

The following dimensionless quantities and correlations were used in this CFD study. 

6.4.2.1 Dimensionless Quantities 


Re r 


co D ^ 

Rem = — (Valensi number) 

4 v 

Mm, max D h , „ , , , , 

(maximum Reynolds number) 


x T 


/ ) j^ e (dimensionless length used for friction-factor plots) 


(6.1) 

(6.2) 

(6.3) 


x* (dimensionless length used for Nusselt number plots) (6.4) 

D/, RePr 


Id = 
If 

Nu_m 

Nu x 


- Ax 


(Darcy friction factor) 


(Fanning friction factor) 

— (mean Nusselt number) 


(6.5) 


( 6 . 6 ) 

(6.7) 


6.4.2.2 Friction Factor and Heat Transfer Correlations 

In order to compare the steady-flow CFD results to the literature, the following correlations were 
selected from Shah and London (1978). The Fanning friction- factor correlation of equation (6.8) is 
attributed to Shah (1978) and applies to laminar, hydrodynamically developing flow, the flow regime for 
the case studied. 


0.674 3.44 'l 

24 -I — , 

J2L + H (6 .g 

(x+) 1/2 1 + 0.000029 (x+) 2 

v 

For steady heat transfer, the correlation of equation (6.9) was selected from Shah and London (1978) 
and is attributed to Stephan (1959). This correlation applies to laminar simultaneously (thermally and 
hydro-dynamically) developing flow which is the flow and thermal regime for the case studied. This 
correlation is valid for a constant wall temperature and Prandtl numbers between 0.1 and 1000. These 
conditions are also true for the cases studied under these CFD steady-state simulations: 
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(6.9) 


Nu m = 7.55 + 


0.024 (x*) 1 14 
1 + 0.0358 (x*)'°- 64 Pr 0 - 17 


This correlation will be referred to in the results section as the “Stephan Nu” or “Stephan (1959)” 
correlation. 

As for the oscillatory-flow cases, the following correlations, of equations (6. 10) and (6. 11), for 
involute-foil friction factor and heat transfer are attributed to Gedeon (see eqs. (3.2) and (3.4) and 
discussion). These correlations were obtained from involute-foil experimental data. The experiments were 
done at Sunpower, Inc., on a NASA/Sunpower oscillating-flow test rig equipped with a microfabricated 
involute-foil regenerator. This regenerator had 42 disks in its stack. The material used for the disks was 
nickel. 


f D =im + o.38Re-° 053 (6.10) 

Re 

This correlation will be referred in the results section as the “Gedeon f Darcy correlation.” Heat transfer 
under oscillatory flow conditions is given by: 

Nu_m = l + 1.97 Pe 0374 (6.11) 


This correlation will be referred in the results section as the "Gedeon Nu_m correlation" 

6.4.2.3 Boundary Conditions 

Steady-state runs: For steady-state runs, the solid temperature was kept constant at 673 K while the 
fluid enters the channel at 660 K. 

Oscillatory-flow runs: The running conditions for the base case examined in the oscillatory-flow 
study are shown in table 6.6. 

TABLE 6.6— BASE CASE FOR OSCILLATORY-FLOW CONDITIONS 


Valensi number, Re,., 0.22885 

Maximum Reynolds number, Re max 49.78 

Frequency, Hz 27.98 

Hydraulic diameter, m 0.000162 

Max mass flux, kg/m 2 -s 6.17215 

Cold end solid B.C Adiabatic 

Hot end solid B.C Adiabatic 

Inlet fluid temperature, cold end, K 293.1 

Inlet fluid temperature, hot end, K 310.2 

Mean pressure, Pa 2500000 

Mean, max velocity, m/s 1 .5488 


6.5 Summary of All Cases Studied 

Tables 6.7 and 6.8 summarize the steady-state and oscillatory-flow conditions, respectively, 
attempted in this study. 
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TABLE 6.7— SUMMARY (STUDY MATRIX) OF STEADY-STATE RUNS 


Geometry 

Material 

Transient or steady 

Re 

3-D 

SS 

steady 

50 

94 

183 

449 

1005 

2213 


TABLE 6.8— SUMMARY (STUDY MATRIX) OF OSCILLATORY FLOW RUNS 


Case Examined 

Geometry 

Material 

Transient or 
steady 

Layer thermal 
contact 
resistance 

Re 

AVt rnax 

Re M 

Comparison with literature, parallel 
plates 

2-D 

Stainless 

steel 

steady 

zero 

Re = 50 


Base case 2-D 

oscillatory 

flow 

zero 

50 

0.229 

Effect of thermal contact 

infinite 

50 

0.229 

Effect of oscillation amplitude 

zero 

150 (3x) 

0.229 

Effect of frequency 

zero 

150 

0.687 (3x) 

Effect of material change 

Nickel 

infinite 

50 

0.229 

Comparison with literature 

3-D 

SS 

steady 

zero 

Re = 50 


Base case 3-D 

oscillatory 

flow 

zero 

50 

0.229 

Effect of thennal contact 

infinite 

50 

0.229 


6.6 CFD Grid-Independence Test and Code Validation 

Three computational domains were identified as good candidates for modeling the problem at hand. 
One was a 2-D domain and the other two were 3-D domains. The questions were: What is a proper grid? 
Will the number of grid locations required to get good results be prohibitively large? From the 
examination of the geometry, one can see that the real involute-foil flow in one channel of one layer 
(fig. 6.4) should approximate that of a parallel-plate geometry. And from the above section we assume 
that the flow is laminar. For this geometry under laminar flow there are theoretical correlations available 
in 2-D. The 2-D computational domain presented above has geometric parameters equivalent to the 3-D 
geometry (same hydraulic diameter, length). It was then natural to select the 2-D domain for studying grid 
independence and code validation 

6.6.1 2-D CFD Grid-Independence-Study Results 

2-D computational domains with four different grid sizes were chosen. The grid sizes (in-the-x- 
directionX in-thc-v-dircction, per layer) are: 20X10, 30X20, 50X20 and 100X40. These are shown in 
figures 6.12 to 6.15 respectively and are summarized in table 6.9. In some cases, the sizes of the grid cells 
in a certain direction are not uniform. A ratio between two successive cell widths is imposed (ratios 
imposed are indicated in the figure labels and in table 6.9), and that allows for denser grids where the 
velocity and temperature gradients are expected to be higher. For channel flow, it is appropriate to have a 
denser grid in the neighborhood of channel walls, and at interfaces between axial-segments where 
discontinuies in wall thermal boundary conditions (or flow areas) occur, since velocity and temperature 
gradients are usually the largest in such regions. 








A 

V 



250 |jm, 20 cells, uniform spacing, 
ratio 1 


Figure 6.12.— The 20X10 2-D grid. 
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81 |jm, 20 cells, ratio 1.1 


250 |_im, 50 cells, ratio 1.1 


Figure 6.14.— The 50X20 2-D grid. 



Figure 6.15.— The 100X40 2-D grid. 


TABLE 6.9— SUMMARY OF GRIDS TESTED IN GRID-INDEPENDENCE STUDY 


Grids/axial-segment 

Number cells 
between plates 

Vertical grid-spacing 
ratio 

Number cells 
along axial 
segment (layer) 

Horizontal or 
axial (segment) 
spacing ratio 

20X10 

10 

1.15 

20 

1, unifonn 

30X20 

20 

1.15 

30 

1 

50X20 

20 

1.1 

50 

1.1 

100X40 

40 

1.1 

100 

1.05 


The above four grid sizes were tested and the results were plotted for friction factor as a function of 
dimensionless length, x + , for Reynolds number, Re = 150 (see fig. 6.16). Similarly, the results were 
plotted for mean Nusselt number, Nu m as a function of the different dimensionless length, x — also for 
Reynolds number, Re = 150 (see fig. 6.17) and for Re = 1000 (see fig. 6.18). The results show poor 
results obtained from the smallest grid (20X10) while little gain in accuracy is achieved by moving from a 
grid size of 50X20 to 100X40 (The 50X20 grid size was eventually selected as the best compromise 
between accuracy and computing resources, to be used following the completion of the grid-independence 
study). 


NASA/CR— 2007-2 15006 


99 





- O - CFD 20x10 — □ — CFD 30x20 — n — CFD 50x20 * CFD 100x40 

Figure 6.16. — Grid independence study: Darcy friction-factors, f_Darcy, as functions of dimensionless length, x + , at 
Reynolds number, Re = 1 50 {for grids/segment of 20X1 0, 30X20, 50X20 and 1 00X40 (horizontal X vertical)}. 



-o- CFD, 20X10 — □ — CFD, 30x20 6 CFD, 50X20 X CFD, 100X40 

Figure 6.17. — Grid independence study: Mean Nusselt numbers, Nu_m, as a functions of dimensionless length, x , at 
Reynolds number, Re = 1 50 {for grids/segment of 20X1 0, 30X20, 50X20 and 1 00X40 (horizontal X vertical)}. 
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x-star 

|- Q - CFD, 20X10 — O — CFD, 30x20 ■ A -CFD, 100X40 | 

Figure 6.18. — Grid independence study: Mean Nusselt numbers, Nu_m, as a functions of dimensionless length, x , at 
Reynolds number, Re = 1 000 {for grids/segment of 20X1 0, 30X20 and 1 00X40 (horizontal X vertical)}. 

6.6.2 CFD Code Validation 

It is best if the CFD solution can be checked, and validated, against experimental data or an exact 
analytical solution of the same problem. However, sometimes it is necessary to do an approximate 
validation if data and analytical solutions are not available. In the case studied, the solution for an 
approximate 2-D grid was chosen as representative of the 3-D channel flow and compared to existing 
theoretical correlations. 

Figures 6.19 and 6.20 show the comparison between the CFD-calculated friction factor and the Shah 
(1978) correlation for Reynolds numbers, Re = 150 and 1000, respectively. At Reynolds number of 150 
the 30X20 and 50X20 friction factors agree well with the correlation. At Reynolds number of 1000 the 
30X20 CFD results seem to be lower than the theoretical correlation at the entrance of the channel but 
agree well downstream where the flow is closer to being fully-developed. Figures 6.21 and 6.22 show the 
comparison between the CFD calculated Nusselt numbers and the Stephan correlation, as reported by 
Shah and London (1978), for Reynolds numbers of 150 and 1000, respectively. Again the agreement is 
generally good except for the entrance region. The entrance region seems to be a difficult area for both 
measurements and computations. Overall there is good agreement with the correlation for both the 
friction factor and the mean Nusselt number and one can expect meaningful results from both the fluid 
flow and the heat transfer analysis from these grids. 
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0 CFD 30x20 — □ — CFD 50x20 - & - Shah (1978) 

Figure 6.19, — Code validation: Darcy friction-factors, f_Darcy, as functions of dimensionless length, x + , at Reynolds 
number, Re = 150; comparisons of CFD computations for grids/segment of 30X20 and 50X20 with Shah (1978) 
correlation. 



— P— CFD 30x20 - a - Shah (1978) 

Figure 6.20. — Code validation: Darcy friction-factors, f_Darcy, as functions of dimensionless length, x + , at Reynolds 
number, Re = 1 000; comparison of CFD computations for grids/segment of 30X20 with Shah (1978) correlation. 
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x-star 


-CFD, 50X20 - A - Correlation, Stephan 


Figure 6.21 . — Code validation: Mean Nusselt numbers, Nu_m, as functions of dimensionless length, x , at Reynolds 
number, Re = 1 50; comparisons of CFD computations for grids/segment of 30X20 and 50X20 with Stephan 
(1959) correlation. 



x-star 


-CFD, 30x20 - * - Correlation, Stephan 


Figure 6.22. — Code validation: Mean Nusselt numbers, Nu_m, as functions of dimensionless length, x , at Reynolds 
number, Re = 1 000; comparisons of CFD computations for grids/segment of 30X20 with Stephan (1 959) 
correlation. 
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6.6.3 Summary: Number and Types of Cells used for Various 2- and 3-D CFD Cases 

Table 6. 10 shows the number and type of cells used in the following study for all 2-D and 3-D cases. 


TABLE 6.10— DESCRIPTION OF THE NUMBER AND TYPE OF CELLS USED. 


Case description 

Zone 

Cell count 

Type of cell 


2-D 6 Layer Grid, 
see figure 6.1 1 

fluid 

6000 

Quadrilateral 

solid 

1800 

Quadrilateral 

total 

7800 


Mesh file size, Mb 

0.67 


Case Description 

Zone 

Cell count 

Type of cell 

3-D straight channel grid, 
see figure 6.10 

fluid 

1404000 

Hexahedral 

solid 

421200 

Hexahedral 

total 

1825200 


Mesh file size, Mb 

372 


Case description 

Zone 

Entry layer cell count 

Repeating unit cell 
count 

Type of cells 

3-D involute-foil channel 
grid, see figure 6.9 

fluid 

582480 

2329920 

Hexahedral 

solid 

93600 

374400 

Hexahedral 

total 

676080 

2704320 


Mesh file size, Mb 

135 

523 



6.7 Results of Two-Dimensional (2-D) CFD Simulations of Involute-Foil Layers 

Table 6.11 shows a summary (study matrix) of all 2-D cases, steady and oscillatory flow, made 
following the grid-independence studies. Steady state was examined first at Reynolds number, Re = 50 in 
order to compare with available correlations. Then oscillatory-flow cases were conducted to examine the 
effects (on friction- factor and Nusselt number) of changing: 1) the thermal contact resistance between the 
6 layers, 2) the oscillation amplitude, 3) the oscillation frequency, and 4) the solid material. 


TABLE 6.1 1— SUMMARY (STUDY MATRIX) OF 2-D CFD SIMULATIONS 


Case examined 

Material 

Transient 
or steady 

Layer thermal-contact 
resistance 

R^max 

Re m 

Comparison with literature, parallel plate 

Stainless 

Steel 

Steady 

zero 

Re = 50 

NA 

Base case 

Oscillatory 

flow 

zero 

50 

0.229 

Effect of thermal contact resistance 

infinite 

50 

0.229 

Effect of oscillation amplitude 

zero 

150 (3x) 

0.229 

Effect of frequency 

zero 

150 

0.687 (3x) 

Effect of material change 

Nickel 

infinite 

50 

0.229 


6.7.1 2-D Steady-State Simulations 

Figures 6.23 and 6.24 show friction factor and Nusselt number comparisons (respectively) for 
Reynolds number, Re = 50. This is also the maximum Reynolds number for the base-case oscillatory-flow 
run. Except for a small region at the entrance, the simulation results agree well with the correlations from 
Shah and London (1978). 
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— Shah (1978) -»-2D Re = 50 

Figure 6.23. — Darcy friction factor, f_Darcy, as functions of dimensionless length, x + : Comparison of values computed 
from 2-D CFD simulations (50X20 grids/segment) with Shah (1978) correlation, at Reynolds number, Re = 50. 



| Stephan (1959) — «— 2D Re = 50 

Figure 6.24. — Mean Nusselt numbers, Nu_m, as functions of dimensionless length, x*: Comparison of values 
computed from 2-D CFD simulations (50X20 grids/segment) with Stephan (1959) correlation, at Reynolds 
number, Re = 50 
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6.7.2 2-D Oscillatory-Flow Simulation 

6.7.2. 1 Base-Case 2-D Oscillatory Flow 

For the oscillatory-flow simulations the base-case forcing-function, at 27.98 Hz, is: 

Mass flux = 6.17215*cos(2*jt*27.98*t + 1.56556) (kg/m 2 -s) (6.12) 

This function is applied at the west (left) fluid boundary. Figure 6.25 shows the variation of the mass flux 
with the crank angle. By monitoring an oscillatory flow variable through several cycles, one notices that it 
takes several cycles until the monitored variable starts varying between the same minimum and maximum 
values. The final condition is called cycle-to-cycle convergence, or is said to have converged to a steady- 
periodic cycle. All the following 2-D oscillatory-flow cases were run until cycle-to-cycle convergence 
was obtained and only after that were the data extracted. Figure 6.26 shows that for the base case it took 
approximately 10 cycles to obtain cycle-to-cycle convergence. The variable monitored in this case is the 
mean fluid temperature in the middle of layer three. 

For characterizing the oscillatory flow, the friction factor is compared with the experimental 
correlation obtained by Gedeon (eq. (6. 10) or (3.2)) for the involute-foil. The friction factor plotted versus 
the crank angle is shown in figure 6.27. The values obtained by the present simulation fall below the 
correlation. This was expected since the correlation was obtained from experimental results on an actual 
involute-foil regenerator while the present 2-D simulation represents an idealized case, with flow through 
a foil-channel that doesn’t flow around foils in adjacent layers (i.e., there are no obstacles in the flow 
path). 



mass flux 

Figure 6.25. — Mass-flux forcing-function as a function of crank angle, CA, degrees, for the base-oscillatory-flow case. 
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Figure 6.26. — Temperature in the middle of layer three as a function of crank angle, CA\ degrees, shows progress 
toward cycle-to-cycle convergence for 2-D base oscillatory-flow case. 



[ ♦ present work * gedeon | 

Figure 6.27. — Darcy friction factors, f_Darcy, as functions of crank angle, CA, degrees; comparisons of values 
calculated from 2-D CFD base-oscillatory-flow case (50X20 grids/segment) with Gedeon involute-foil correlation 
of equation (6.10) and (3.2). 
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In order to characterize the heat transfer that takes place during an oscillatory-flow run, the mean 
Nusselt number is plotted with respect to the crank angle, in figure 6.28. For comparison the experimental 
correlation, also obtained by Gedeon (eq. (6. 1 1) or (3.4)), for the mean Nusselt number is used. However 
the correlation obtained by Gedeon represents a mean over the length of a stack of 42 layers tested at 
Sunpower Inc. The present work focuses only on the region from the middle of layer 3 to the middle of 
layer 4. So the length over which the Nusselt number is averaged in the present work is equal to the 
thickness of one layer only. That is done in order to stay away from the ends where entrance effects can 
distort the results. On the other hand, calculating a mean Nusselt number over the length of 6 layers 
would not have been representative of the actual geometry. Furthermore, experimental testing done at the 
University of Minnesota looked at the mean Nusselt number calculated between two layers similar to 
what has been done in the present work. The shape of the mean-Nusselt-number, or Nu_m, plot versus the 
crank angle is different from the Gedeon correlation and this difference arises from how the Nusselt 
number is averaged. The experiments done at the University of Minnesota (UMN) show a Nusselt number 
curve similar to the present work. However the comparison with Gedeon correlation is useful for the 
maximum Reynolds number regions that are located around 90° and 270° crank angle, where flow rates 
are maximum in the two directions. At these locations the 2-D analysis lies slightly below the correlation. 

By integrating the fluid enthalpy crossing a plane section in the middle of layer three over the whole 
cycle, one can calculate the net enthalpy loss over the cycle. If one integrates the solid-conduction heat 
transfer at the middle of layer three over the whole cycle, a net-conduction loss over the cycle can be 
obtained. Since both losses are crossing the plane at the middle of layer three they can be added together 
to obtain a total axial heat loss over the cycle. For the base-oscillatory-flow case, table 6.12 shows these 
2-D CFD heat-loss results. 



crank angle 


present work Gedeon | 

Figure 6.28. — Mean Nusselt numbers, Nu_m, as functions of crank angle, CA, degrees; comparison of values 
calculated from 2-D CFD base-oscillatory-flow case (50X20 grids/segment) with Gedeon involute-foil correlation 
of equation (6.1 1 ) or (3.4) {CFD assumes perfect thermal contact between layers, or zero thermal contact 
resistance (TCR)} 
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TABLE 6.12— 2-D CFD BASE-OSCILLATORY-FLOW-CASE ENTHALPY, 
CONDUCTION AND TOTAL-AXIAL HEAT LOSSES 



Enthalpy loss, 
W 

Conduction loss, 
W 

Total loss, 
W 

Base case 

1.722 

1.174 

2.896 


6.7.2.2 Effect of Changing the Thermal Contact Resistance 

In order to study the effect of thermal contact resistance (TCR) between the layers a change was made 
from the zero TCR of the base case to an infinite TCR condition at the interfaces between the layers (or 
from perfect thermal contact to perfect thermal insulation, between solid layers). The effect of thermal 
contact resistance is important because in reality the contact between the layers is not perfect. This added 
contact resistance impedes the solid conduction from layer to layer. This contact resistance between the 
layers causes discontinuities in the solid-wall temperature profile between the hot and cold sides of the 
interface. The changed wall-temperature profile should affect the heat transfer between the wall and the 
fluid which, in turn, should change the plot of the Nusselt number. However the friction factor is not 
affected (as expected). Figure 6.29 shows a comparison of the Nusselt number behavior between the base- 
case, zero-TCR and the infinite-TCR cases. 

The infinite TCR (adiabatic contact) has caused the Nusselt number to rise, especially in the regions 
of low Reynolds numbers close to where the flow reverses (near 180° and 360° crank angle). This was 
expected. When no TCR is present, the solid wall temperature on each side of the contact between the 
solid layers is the same. When the infinite TCR is introduced, a temperature difference in the solid wall 
between the two sides of the contact results. The fluid flowing in the channel past the contact between the 
two solid layers sees the discontinuity in the wall temperature profile. The increased delta-T between the 
wall and the fluid causes an increase in the heat transfer and that is reflected in the higher Nusselt number. 
At lower Reynolds numbers, the fluid has more time to absorb the heat and the effect of the infinite 
thermal contact resistance is more pronounced. However, when the flow stops, as it does when it switches 
direction, the temperature between the fluid and the solid equalizes and delta-T becomes very small, 
tending to zero. That introduces a discontinuity in the calculation of the Nusselt number and figure 6.29 
shows that discontinuity at 180° and 360° of crank angle. The other change that has happened is related to 
the difference between the cooling and the heating parts of the cycle. The cooling half happens from zero 
to 1 80° of crank angle when the flow goes from the cold side to the hot side, and the heating half happens 
from 180° to 360°. When zero thermal contact resistance was present (base-case) the mean-Nusselt- 
number curves for the two halves looked the same. However with infinite thermal contact introduced, the 
heating half of the cycle shows a higher mean Nusselt number. 

In terms of heat loss, table 6. 13 shows the results of integrating the enthalpy loss and the conduction 
loss over the whole cycle. In the case of the solid conduction, increasing the TCR from zero (base case) to 
infinite has resulted in a reduction of 54.7%. The enthalpy loss has increased by 13.8%. However the total 
loss has decreased by 14.0%. This suggests that is a good idea to increase the thermal contact resistance in 
order to reduce the heat loss. 
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Figure 6.29. — Mean Nusselt numbers, Nu_m, as functions of crank angle, CA, degrees; comparison of values 

calculated from 2-D CFD oscillatory-flow cases (50X20 grids/segment) for (1) zero TCR (perfect thermal contact, 
base case) and (2) infinite TCR (adiabatic contact or perfect insulation). 



TABLE 6.13— HEAT-LOSS COMPARISON OF ZERO-TCR BASE CASE TO INFINITE-TCR CASE 



Enthalpy loss, 
W 

Change 

Conduction loss, 
W 

Change 

Total loss, 
W 

Change 

Base Case 

1.722 


1.174 


2.896 


Infinite TCR 

1.960 

13.8% 

0.531 

-54.7% 

2.491 

-14.0% 


6.7.2.3 Effect of Changing the Oscillation Amplitude 

The increase of the oscillation amplitude results in an increase of the maximum Reynolds numbers at 
90° and 270° CA. Figure 6.30 shows how the behavior of the friction factor has changed due to an 
increase in oscillation amplitude by a factor of three. The friction factor has dropped and that is consistent 
with the steady-state simulations. When the Reynolds number is increased, the friction factor becomes 
smaller. Figure 6.3 1 shows how the mean Nusselt number changes when the oscillation amplitude (and, 
consequently, the maximum Reynolds number) are increased by a factor of three. The figure shows that 
Nusselt number stays much the same. When the hot or the cold fluid is pushed deeper into the stack the 
temperature difference between the wall and the fluid should increase. That represents the driving 
potential for the heat transfer. The figure tells us that the heat flux has also increased in such a way that 
the heat transfer coefficient has stayed about the same. The increased flow due to the increase in 
amplitude of the mass-flow would be responsible for sustaining the increased wall heat flux. 

In terms of heat loss, table 6. 14 shows the results of integrating the enthalpy loss and the conduction 
loss over the whole cycle, for the base case and the increased-oscillation-amplitude case. 
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TABLE 6.14— HEAT LOSS COMPARISON OF BASE CASE (Re max = 50) 
AND INCREASED-OSCILLATION-AMPLITUDE CASE (Re max = 150) 



Enthalpy loss, 
W 

Change 

Conduction loss, 
W 

Change 

Total loss, 
W 

Change 

Base Case 

1.722 


1.174 


2.896 


Remax = 150 

18.249 

10. 6x 

0.956 

-18.6% 

19.204 

6.6x 


A reduction in conduction loss occurred after the oscillation amplitude was increased by a factor of 
three. However, due to the higher flow, the enthalpy loss increased by a factor of 10.6. The effect on the 
total heat loss is an increase by a factor of 6.6. The reduction in the axial conduction is attributed to the 
higher heat flow from gas to metal due to the higher instantaneous gas mass flow. This resulted in heating 
the solid faster in the axial direction and thus lower (instantaneous) axial heat conduction resulted. 



CA 

— 0— Re_max =150 — □ — Re max - 5Q 


Figure 6.30. — Darcy friction factors, f_Darcy, as functions of crank angle, CA, degrees; comparisons between 2-D 
CFD oscillatory-flow cases (50X20 grids/segment): base-case (Re max = 50) and increased-oscillation-amplitude 
case (Re max = 150). 
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Figure 6.31 . — Mean Nusselt numbers, Nu_m, as functions of crank angle, CA, degrees; comparisons between 2-D 
CFD oscillatory-flow cases (50X20 grids/segment): base case (Re max = 50) and increased-oscillation-amplitude 
case (Re max = 150). 


6.7.2.4 Effect of Changing the Oscillation Frequency 

The effect of increasing the oscillation frequency was also studied. As shown in the CFD oscillatory- 
flow test matrix (table 6.8, summary of oscillatory-flow runs), a frequency that is three times the base- 
case frequency was chosen for testing this effect. Therefore, frequency increased from 27.98 to 84Hz 
while the Valensi number, Re m or Va, increased from 0.229 to 0.687. A change in frequency would alter 
both the fluid flow and the heat transfer behavior. Figure 6.32 shows the behavior of the friction factor 
when the oscillation frequency is increased three times. It can be seen that at low Reynolds numbers, in 
the vicinity of where the flow reversal takes place (-180° and 360°), the friction factor deviates from the 
one at the base-case frequency. When the fluid is decelerating prior to 180°, the friction factor becomes 
smaller. At crank angles immediately after flow reversal, the friction factor becomes larger. Similarly, the 
shape of the Nusselt number curve in figure 6.33 changes when the frequency is increased. To 
summarize, when the fluid is accelerating, the Nusselt number is higher than the base-case number and 
when the fluid is decelerating, the Nusselt number becomes lower than the base-case number. 

In terms of heat loss, table 6. 15 shows the results of integrating the enthalpy loss and the conduction 
loss over the whole cycle for the base-case (Re max = 50, Re m = 0.229), the increased-oscillation-amplitude 
case (Re max = 150, Rc„, = 0.229), and the increased-oscillation- frequency case (Re max = 150, Rc m = 
0.687). 
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TABLE 6.15— HEAT LOSS COMPARISON OF BASE CASE, (Re max = 50, Re m = 0.229), THE INCREASED-OSCILLATION- 
AMPLITUDE CASE (Re max = 150, Re,,, = 0.229) AND THE INCREASED-OSCILLATION-FREQUENCY 
CASE (Re max = 150, Re m = 0.687) 



Enthalpy loss, 
W 

Change 

Conduction loss, 
W 

Change 

Total loss, 
W 

Change 

Base case 

Re m ax = 50, Re l0 = 0.229 

1.722 


1.174 


2.896 


Inc. osc. amplitude 
Re m ax = 150, Re m = 0.229 

18.249 

10. 6x 

0.956 

-18.6% 

19.204 

6.6x 

Inc. osc. frequency 
Re max = 150, Re m = 0.687 

13.466 

7.8x 

1.009 

-14.1% 

14.474 

5. Ox 



|^^Re_max =150, Re_co = 0.229 * Rejnax =150, Re_u) = 0,687 | 

Figure 6.32. — Friction-factor comparison between the base case (Re max = 50, Re w = 0.229) and the increased- 
oscillation-frequency case (Re max = 150, Re u = 0.687). 
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Figure 6.33. — Mean Nusselt number comparison between the base case (Re max = 50, Re w = 0.229) and the 
increased-oscillation-frequency case (Re ma x = 150, Re^ = 0.687). 

A 14% reduction in the axial conduction took place while an increase (by a factor of 7.8) in the 
enthalpy flux occurred when compared to the base case. The net effect is about a factor of 5.0 increase in 
the regenerator axial heat loss. Again, the reduction in axial conduction is attributed to the higher heat 
flow from gas to metal due to the higher instantaneous gas mass flow. This resulted in heating the solid 
faster in the axial direction and, thus, lower (instantaneous) axial heat conduction resulted. 

6.7. 2.5 Effect of Changing the Solid Material 

The impact of solid-material properties on performance are of interest. Pure metals are known to have 
higher conductivity than alloys. Nickel was chosen for fabrication of the prototype test regenerator, due to 
limitations of cost and time — since a nickel regenerator could be fabricated via LiGA alone. However, 
nickel has a higher conductivity than stainless steel (the preferred material) by about 5.5 times (see 
table 3.10). For this comparison of nickel and stainless steel materials, the contact between layers has 
been set to infinite thermal contact resistance (TCR); thus the reference stainless-steel case for this 
material study was not the base case (which had zero TCR). Figure 6.34 compares the friction-factor 
behavior and, as expected, no change in the friction factor has been detected upon changing the material 
from stainless steel to nickel. Figure 6.35 shows the behavior of the mean Nusselt number upon changing 
the material. Because of the higher conductivity of the nickel the wall temperature profile was flatter than 
that for stainless steel for the length of one layer, between two interfaces of infinite thermal contact 
resistance. That should cause a change in the heat transfer. Note that figure 6.35 compares the cases of 
nickel with infinite TCR (adiabatic contact), stainless steel with infinite TCR (adiabatic contact) and the 
case of stainless steel with zero TCR (perfect contact, the base case). 
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—□—Nickel * Stainless Steel 


Figure 6.34. — Friction-factor comparison between nickel and stainless steel materials, both with infinite TCR between 
layers (adiabatic contact). 



—^"Nickel, adiabatic contact —♦—Stainless Steel, adiabatic contact * ■ Stainless steel, perfect contact | 

Figure 6.35. — Mean Nusselt number comparison between nickel and stainless steel, both with infinite TCR (adiabatic 
contact) between layers, and the base case (stainless steel, zero TCR or perfect contact). 
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The results show that, for infinite TCR between sections, the mean Nusselt number has increased 
overall for nickel compared to that for stainless steel, especially at low Reynolds numbers (i.e., near flow 
reversals at -180° and 360°). Again, this is explained by the fact that within a given layer the temperature 
profile is flatter resulting in higher temperature differences between the wall and the fluid. The lower 
conductivity material can maintain a steeper temperature profile within one layer. This steeper profile is 
closer to the bulk-temperature profile of the fluid, leading to smaller temperature differences between 
fluid and solid. A comparison of the axial heat losses for nickel and stainless steel is given in table 6. 16 
(both with infinite TCR or adiabatic contact at the interfaces between layers). 


TABLE 6.16— HEAT LOSS COMPARISON BETWEEN STAINLESS-STEEL AND NICKEL MATERIALS, 
BOTH WITH INFINITE TCR (OR ADIABATIC CONTACT) AT INTERFACES BETWEEN LAYERS (OR DISKS) 



Enthalpy loss, 
W 

Change 

Conduction loss, 
W 

Change 

Total loss, 
W 

Change 

Infinite TCR and SS 

1.960 


0.531 


2.491 


Infinite TCR and Nickel 

1.862 

-5.1% 

0.724 

+36.3% 

2.586 

+3.8% 


In this case a 36.3% increase in the axial conduction took place while a decrease (by 5.1%) in the 
enthalpy flux occurred when nickel is compared with stainless steel, both cases with infinite TCR 
between layers. The net effect is about a 3.8% increase in the regenerator heat loss using nickel. The 
increase in the axial conduction is attributed to the higher thermal conductivity of nickel (compared to 
stainless steel). 


6.8 Results of 3-D CFD Straight-Channel-Layers 
Approximation of Involute-Foil Layers 

The focus of this section is the presentation and discussion of the results for the 3-D straight-channel- 
layers (fig. 6.10) simulations. As mentioned earlier, two types of 3-D simulations were developed. One is 
the 3-D straight-channel (an approximation of an involute-foil) and the other is the 3-D involute-foil 
channel. The grid for the involute-foil problem is quite dense and it was not feasible for use in oscillatory- 
flow simulations. The 3-D straight channel grid was subjected to oscillatory flow boundary conditions 
and the results are to follow. 

As with the actual microfabricated involute-foil geometry, the flow channels alternate in orientation 
from layer to layer. Figure 6.10 shows that at the entrance, the channels are horizontal for this view and 
on the second layer, the channels are vertical. This results in a small contact area between layers. 
Flowever, the walls of the second layer disturb the flow leaving the first layer by forcing the flow to 
reorient itself to the new channels. That should result in an increase of friction factor. Flowever the heat 
transfer should improve because of the boundary-layer disturbance. The question of how the thermal 
contact between the layers affects the heat transfer is again important. Allowing maximum solid 
conduction (zero TCR) through the interfaces between layers again forms the base case for these straight- 
channel simulations. The effect of setting an infinite thermal contact resistance between the layers is then 
studied and compared with the base case and also with the 2-D simulation. Table 6.17 shows the 
conditions, or study matrix, for the 3-D simulations performed using the straight-channel geometry. 


TABLE 6.17— SUMMARY, OR STUDY MATRIX, FOR 3-D STRAIGHT-CHANNEL-LAYERS RUNS 


Case examined 

Material 

Transient or 
steady 

Layer thennal-contact 
resistance (TCR) 

Re 

lvt max 

Re w 

Comparison with literature 

SS 

steady 

Zero 

50 

0.0 

Base case 

oscillatory flow 

Zero 

50 

0.229 

Effect of thermal contact resistance 

oscillatory flow 

Infinite 

50 

0.229 
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6.8.1 Steady-Flow Simulation (3-D Straight-Channel Layers) 

A 3-D steady-state simulation at a Reynolds number of 50 was performed in order to compare with 
the literature and with the 2-D results. In order to compare with the 2-D results, and with the literature, the 
Darcy friction factor is plotted against the hydrodynamic axial coordinate, x + , which is dimensionless. 

The same correlation from Shah (1978), discussed earlier, is used. In the 2-D simulation, the Shah 
correlation was below the 2-D CFD results in the entry section and then matched well with the results. 
Figure 6.36 compares these 3-D steady-flow results with the 2-D results and the correlation from Shah. 

As can be seen in figure 6.36, the 3-D results agree well with the 2-D results for the first layer. Upon 
entering the second layer, the flow encounters some resistance from the perpendicularly-oriented second 
layer. That is where the friction factor goes up and departs from the agreement with the 2-D results. Then 
the flow tries to settle again until it encounters another geometry change upon entry into the third layer. 

As it moves through the stack of layers, the behavior of the fluid flow settles into periodicity, with small 
increases in friction factor upon entering each layer and with an average value above the 2-D prediction. 
This behavior was expected and the simulation provided an answer regarding the magnitude of the 
friction factor increase. 

In order to characterize the heat transfer, the mean Nusselt number is plotted with respect to the 
dimensionless thermal axial coordinate, x , in figure 6.37 and is compared to the results from the 2-D 
simulations and the correlation developed by Stephan (1959). As discussed earlier, the alternating 
orientation of the layers is expected to improve the heat transfer, relative to 2-D and uniform- channel 
flow. That should result in higher 3-D Nusselt number values at each flow-channel discontinuity; this can 
be observed in figure 6.37 . 

As can be seen in figure 6.37, the 3-D results agree well with the Stephan (1959) correlation and the 
2-D simulation for the first layer. However, the Nusselt number increases for the following layers. This is 
expected because of the disturbance in the thermal boundary layer introduced by the changing orientation 



X+ 

—♦—Shah (1978) -«-3D Re = 50 — *-2D Re = 50 

Figure 6.36. — 3-D straight-channel-layers, steady-flow, friction-factor comparisons with 2-D results and Shah (1978) 
literature correlation — all at Reynolds number, Re = 50. 
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—♦—Stephan (1959) — 3D Re = 50 -*-2D Re = 50 

Figure 6.37. — 3-D straight-channel-layers, steady-flow, mean-Nusselt-number comparisons with 2-D results and 
Stephan (1959) literature correlation — all at Reynolds number, Re = 50. 


of the channels. After several layers, the heat transfer settles into a periodic behavior with an average 
Nusselt number higher than the 2-D simulation. The small irregularities in the 3-D and 2-D results 
represent interpolation error due to the location of the data extraction planes with respect to the grid. The 
CFD software stores data at the centers of cells. If data are requested at a location between the cell 
centers, the software performs a linear interpolation, which introduces a small error. 

6.8.2 Oscillatory-Flow Simulation (3-D Straight-Channel Layers) 

6.8.2. 1 Base-Case 3-D Oscillatory Flow 

For the 3-D straight-channel oscillatory-flow simulations, the same forcing-function for the mass flux 
as in the 2-D simulations was used (see eq. (6. 12)). As in the 2-D simulations, it took about 10 cycles to 
establish cycle-to-cycle convergence. The expectation is again that both the friction factor and the Nusselt 
number would be higher. 

As for the 2-D simulations, the friction factor is plotted (see fig. 6.38) against the crank angle in 
degrees and compared to the 2-D case and the experimental correlation for involute-foil by Gedeon 
(eq. (6.10) or (3.2)). 

As expected, the 3-D friction-factor values are higher than the 2-D results at all crank angles and are 
more in line with the experimental involute-foil correlation values from Gedeon (eq. (6.10) or (3.2)). The 
mean-Nusselt-number results determined from the 3-D simulation are compared (see fig. 6.39) with the 
2-D results and the experimental involute-foil correlation from Gedeon (eq. (6.11) or (3.4)). As mentioned 
earlier, the correlation from Gedeon averages the Nusselt number over the length of a stack of layers 
while the present work uses only the length of one layer to obtain the mean Nusselt number. As expected, 
the Nusselt number values are higher than the values for the 2-D parallel-plate simulation and also higher 
than the Gedeon correlation given by equation (6.1 1) or (3.4). 
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—♦— 2D base case Gedeon m 3D base case 

Figure 6.38. — 3-D straight-channel-layers, oscillatory-flow, friction-factor comparison to 2-D base-case oscillatory-flow 
results and the Gedeon (eg. 6.10) experimental involute-foil correlation. 



—♦—2D base case ^ — Gedeon * 3D base case 

Figure 6.39. — 3-D straight-channel-layers, oscillatory-flow, mean-Nusselt-number comparison to 2-D base-case 
oscillatory-flow results and the Gedeon (eq. (6.1 1 ) or (3.4)) experimental involute-foil correlation. 
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6.8.2.2 Effect of Changing the Thermal Contact Resistance (3-D Straight-Channel-Layers) 

As in the 2-D work, in order to study the effect of thermal contact resistance between the layers, an 
infinite thermal-contact-resistance (TCR) condition was imposed at the interfaces between the layers 
(replacing the zero-TCR of the base case). The expectation is that the friction factor will not change when 
compared with the perfect contact case. However, the mean Nusselt number should behave similar to how 
it behaved when the same condition was imposed in the 2-D study. That is, the Nusselt number should 
increase overall, especially at the low Reynolds numbers that are encountered when the flow switches 
direction (near 180° and 360°). Figure 6.40 shows a comparison of the mean-Nusselt-number behavior 
among the 3-D zero-TCR (perfect contact) and the infmite-TCR (adiabatic contact) cases and the infmite- 
TCR (adiabatic contact) 2-D case. When compared to the 3-D zero-TCR case (i.e., perfect contact), the 
infmite-TCR (adiabatic contact) case has higher mean Nusselt number, especially in the regions of lower 
Reynolds numbers close to where the flow reverses. This was expected and it is similar to what happened 
in 2-D when the infinite thermal contact resistance condition was imposed. 



2D adiabatic contact ■ 3D adiabatic contact ♦ 3D perfect contact 

Figure 6.40. — Mean-Nusselt-number comparison for 3-D straight-channel-layers, oscillatory-flow cases, with zero- 
TCR (perfect contact) and infinite-TCR (adiabatic contact) between layers — and for the 2-D oscillatory-flow 
infinite-TCR (adiabatic contact) case. 
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6.9 Results, 3-D Steady-Flow Simulation of Involute-Foil Layers 

In this section, the results for the 3-D involute-foil-layer simulations are presented and discussed. 

This 3-D computational domain resembles more closely the actual microfabricated design in terms of the 
shape of the channel. However, the attempt to capture the geometry of the channel better resulted in a 
dense grid. Even after reducing the length of the stack to only two layers, the cell count was close to 
2.7 million. As explained earlier, one can use a two-layer-long unit repeatedly by taking the velocity and 
temperature profile from the outlet and applying it back to the inlet in order to simulate a longer stack. 
The drawback of this technique is that oscillatory-flow simulation is not possible. Instead of one 
oscillatory-flow run, simulations have been performed using steady-state conditions at several Reynolds 
numbers. 

6.9.1 Summary of 3-D Involute-Foil Runs 

The 3-D involute-foil grid actually consists of two types of grids, one half layer for the entry and a 
repeating unit consisting of a one-half-layer entry, a full layer, and a one-half-layer exit (see fig. 6.9). 
Cutting the first layer in half allows for the two grids to line up properly when passing the boundary 
conditions (profile) from one grid to the other. By the same token, two repeating units will line up 
properly so that the boundary profile can be passed, where the geometry is continuous. As shown in 
figure 6.9, the grid captures one full involute-foil channel in the middle. The other channels are also 
simulated as full because of the periodic boundary conditions applied to the sides. In fact, when 
periodicity is considered, this grid simulates a full ring of channels. The contact between the layers is 
perfect; the thermal contact resistance is zero. 

Meshing this geometry presented a considerable challenge. Even though the channels were insured to 
line up properly, in order not to get discontinuities when passing the boundary profile, one has to ensure 
that the cell coordinates of the exit face match the cell coordinates of the inlet face. The ends of the 
channel with their round shape did not allow for a structured mesh and an unstructured mesh had to be 
employed. However, the generation of the unstructured mesh is harder to control when it comes to 
matching cell coordinates between inlet and outlet faces. This required linking the inlet and the outlet 
faces, which is a tedious process. Table 6.18 shows a summary of simulations performed using the 3-D 
involute-foil grid. 


TABLE 6.18.— SUMMARY, OR STUDY MATRIX, OF 3-D 
STEADY-FLOW INVOLUTE-FOIL SIMULATIONS 


Geometry 

Material 

Transient or steady 

Layer thermal contact 

Re no. 

3-D Involute- 
foil channel 

SS 

Steady 

Zero 

50 

94 

183 

449 

1005 

2213 


When the 2-D steady-flow simulations were performed, it was necessary to impose boundary 
conditions, as summarized in table 6.19. 


TABLE 6.19— INLET BOUNDARY CONDITIONS FOR THE 
SIMULATED REYNOLDS NUMBERS, RE. OF TABLE 6.18 


Case 

B.C. in 
CFD 

Inlet temp, 
K 

Wall temp, 
K 

Velocity inlet, 
m/s 

Re no. 

1 




4.7 

50 

2 




8.9 

94 

3 

Velocity-inlet 

660 

673 

17.3 

183 

4 

42.4 

449 

5 




94.9 

1005 

6 




209.0 

2213 
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6.9.2 Steady-State, Re = 50, Friction-Factor and Nusselt Number Comparisons 

As for the 2-D parallel-plate and the 3-D straight-channel-layers, a 3-D involute-foil-layers simulation 
at Reynolds number 50 was performed. In order to compare with the 2-D results and with the literature, 
the Darcy friction factor is plotted (see fig. 6.41) against the dimensionless hydrodynamic axial 
coordinate, x + . The same literature correlation, Shah (1978), is used. 

The 3-D involute-foil-layers simulation shows a variation in friction factor (the saw shape) similar to 
the 3-D straight-channel layers, as expected. One thing to keep in mind is that the length of the involute- 
foil layer is 15 pm longer in the flow direction than the 3-D straight-channel layer. While work was in 
progress on this project it was learned that the actual fabricated layers were shorter than originally 
intended. The 3-D straight channel was adapted to the shorter length and simulations were performed that 
way. However the 3-D involute-foil layer length simulated was kept at the original length. The above 
comparison captures this difference graphically by showing that the layer-to-layer rise in friction factor 
for the 3-D involute-foil layers happens after the rise shown by the 3-D straight-channel layers. 

In order to characterize the heat transfer, the 3-D involute-foil mean Nusselt number is plotted in 
figure 6.42 with respect to the thermal axial coordinate and is compared to the results from the 2-D 
simulation, 3-D straight-channel simulation, and the correlation developed by Stephan. The variation in 
mean Nusselt number, Nu_m, is similar to that encountered for the 3-D straight channel grid. The out of 
step layer to layer transition is again explained by the difference in the length of the layer. 



'Shah (1978) straight channel ^^^™3D involute channel 


Figure 6.41 . — Steady-state friction-factor comparison at Reynolds no., Re = 50 as a function of dimensionless length . 
Calculated via the 3-D involute-foil and straight-channel layer simulations, the 2-D parallel-plate simulation — and 
compared with the Shah (1978) correlation. 
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X* 

♦ Stephan (1959) — A — 2D — X— 3D straight channel — o— 3D involute channel 

Figure 6.42. — Steady-state mean-Nusselt-number comparisons at Reynolds number, Re = 50 as a function of 

dimensionless length. Calculated via the 3-D involute-foil and straight-channel layer simulations, the 2-D parallel- 
plate simulation — and compared with the Stephan (1959) correlation. 


6.9.3 Summary: Steady 3-D Involute-Foil Friction-Factor Results for All Reynolds Numbers 

Simulations at other values of Reynolds numbers have been performed in order to determine the 
variation of the friction factor with the Reynolds number. A total of six Reynolds numbers have been 
attempted. The results are presented in groups of three from lowest to highest. 

Figure 6.43 compares the variation of the friction factor at Reynolds numbers 50, 94, and 183. The 
friction factor is lower at higher Reynolds numbers, as expected. As in previous simulations, the regular 
friction factor increases observed at the transition from layer to layer are present at the other Reynolds 
numbers. The fact that the transitions don’t line up is an artifact of the plotting versus the dimensionless 
axial coordinate, x+, which includes the Reynolds number in the denominator. Figure 6.44 presents the 
variation of the friction factor versus the actual axial coordinate in meters. One can see that the transitions 
between the layers line up now. However such a representation in terms of the actual dimension is not 
useful when a comparison with other work or an experiment is to be done. The results for the next three 
Reynolds numbers are presented in figure 6.45. Similar to the previous three Reynolds number 
simulations, the friction factor decreases as the Reynolds number increases. Another important 
observation is that as the Reynolds number increases, it takes more layers for the friction-factor variation 
to flatten out. In the case of Reynolds number = 2213 the simulation has been run up to 14 layers while 
for Reynolds number = 50 only six layers appear to be enough for the friction-factor trend to flatten. 
Averages of the friction factor have been taken from the last layers and plotted against the Reynolds 
number and compared with results from experiments and correlations in figure 6.46. 
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Figure 6.43. — Steady 3-D involute-foil simulation friction-factor at different Reynolds numbers, Re = 50, 94, and 
183 — as a function of the dimensionless length, x + . 



x(m) 


| — Re = 50 — Re = 94 —A— Re = 183 | 

Figure 6.44. — Steady 3-D involute-foil simulation friction-factors at different Reynolds numbers, Re = 50, 94, and 
183 — as a function of the actual length, x. 
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Figure 6.45. — Steady 3-D Involute-foil simulation friction factors at different Reynolds numbers, Re = 449, 1005, and 
2213 — as a function of the dimensionless length, x + . 




Figure 6.46. — 3-D involute-foil-simulation friction factor (CFD-CSU-8) 
compared with UMN large-scale measurements (8 and 10 layers), 
Gedeon involute-foil experimental correlation (GA-42-1, eq. (3.1), 
original stacking) and Kays and London (1964) staggered-plate 
correlation. 


NAS A/CR — 2007 -2 1 5006 


125 


In figure 6.46, the results obtained from the 3-D involute-foil simulations have been compared with 
experiments done by the University of Minnesota (UMN) on eight- and ten-layer large-scale stacks (see 
section 5.3.5). The same graph contains a correlation developed by Kays and London (1964) for 
staggered-plate heat exchangers and a correlation developed by Gedeon (eq. (3.1), for original stacking) 
based on involute-foil experiments performed in the NASA/Sunpower oscillating-flow test rig. The 
results from the current simulations match very well with the UMN experiments. Both the CFD and UMN 
data match very well with equation (3.1) at the low end of the Reynolds number range (—100 to 200). 
However equation (3.1) provides higher f values at the high end of the Reynolds number range (-1000). 
This is attributed to the roughness associated with the EDM cutting process, shown earlier in figure 3.7. 

6.9.4 Summary: Steady 3-D Involute-Foil Mean-Nusselt-Number Results for All Reynolds 
Numbers 

The steady 3-D involute-foil simulations performed at the six Reynolds numbers were also analyzed 
for heat transfer. The mean Nusselt number was plotted versus the dimensionless thermal axial 
coordinate, x , which also includes the Reynolds number in the denominator. This causes the plots to 
appear compressed in the axial direction, as Reynolds number increases. Figure 6.47 shows the mean 
Nusselt number variation for the first three Reynolds numbers, 50, 94, and 183. The mean Nusselt 
number increases as the Reynolds number increases. There is a good match for the first layer where the 
graphs appear to overlap. However, for the following layers, the mean Nusselt number, Nu_m, settles at 
higher values for increased Reynolds numbers. The plots appear to flatten as they advance axially and that 
can be associated with the flow reaching a more thermally-developed condition. As with the friction 
factor, at higher Reynolds numbers it takes more layers for the plot to flatten. Figure 6.48 shows the 
results for the next three Reynolds numbers, 449, 1005, and 2213. The behavior of the mean Nusselt 
number follows the same behavior of being the same in the first layer but increasing as the Reynolds 
number increases. 



X* 

pi- Re = 50 — Re = 94 — *— Re = 183 

Figure 6.47. — Steady 3-D involute-foil simulation mean-Nusselt-numbers at different Reynolds numbers, Re = 50, 94, 
and 183 — as functions of dimensionless length, x*. 
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Figure 6.48. — Steady 3-D involute-foil-simulation mean Nusselt numbers at different Reynolds numbers, Re = 449, 
1005, and 2213 — as function of dimensionless length, x* 


6.10 CFD Simulation Conclusions 

This section provides the conclusions to the work described in the previous CFD sections (Starting 
with section 6.0 “Analysis Tools and CFD Results”). A multitude of simulations were run using three 
main geometric representations of involute foils: 2-D parallel-plate, 3-D straight-channel layers, and 3-D 
involute-foil layers. Both steady-state and oscillatory-flow simulations were carried out for the 2-D and 3- 
D straight-channel geometries. Only steady-state simulations were considered practical for the 3-D 
involute-foil geometry (because of grid size and CPU time required). The results focused mainly on two 
quantities of interest, the friction factor and the mean Nusselt number. These quantities were compared 
for various simulation cases, experimental correlations from the literature, and experimental data obtained 
at UMN (large-scale tests) and Sunpower Inc (actual-scale tests). The simulation results were based upon 
the whole flow length for the steady-state cases. In the case of the oscillatory-flow simulations, where the 
computational domain has only six layers, the quantities of interest were averaged only for the length of 
one layer from the middle of layer three to the middle of layer four. 

6.10.1 Conclusions: 2-D Parallel-Plate Simulations (Steady and Oscillating Flow) 

The 2-D parallel-plate steady-state simulations showed good agreement with the correlations from the 
literature. That provided confirmation that the chosen grid refinement is able to resolve the fluid flow and 
the heat transfer for the laminar flow regime. 

A base-line case was chosen for this study using helium as the working fluid at pressure = 2.5E06 Pa 
(25 atm), and stainless steel solid material with zero thermal contact resistance (TCR) between the 6 
layers. The oscillation frequency and amplitude chosen resulted in maximum Reynolds number, 

Re ma x = 50, and Valensi number, Rc m = 0.229. For this base case the enthalpy loss computed throughout 
the cycle = 1.722 W while the axial conduction loss = 1.174 W with a total of 2.896 W. This total value 
must be minimized for optimum regenerator design. 
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Changing the solid thermal contact resistance from zero to infinity between the layers had no effect 
on the friction factor, as expected. However, the mean Nusselt number increased overall with the increase 
more pronounced at the crank angles where the flow switches direction. The infinite thermal contact 
resistance between the layers blocked solid conduction from layer to layer. That resulted in a higher 
temperature difference between the solid and the fluid with an increase in the heat transfer, as represented 
by the higher Nusselt numbers. Comparing the results of this case with the base case, greater than 54.7% 
reduction in the axial conduction took place (as expected) while an increase (by 13.8%) in the enthalpy 
loss occurred. The net effect was about a 14% reduction in the regenerator axial heat loss. This confirms 
the finding of the earlier 1-D analyses shown in section 3. 6. 5. 3., i.e., increasing the TCR between the 
layers will result in a reduction of the axial conduction accompanied with an increase in the enthalpy heat 
loss. 

Increasing the velocity amplitude of the oscillation (by a factor of 3) resulted in lower friction factors. 
The increased amplitude means higher maximum Reynolds number (=150) and lower friction factors due 
to the fact that the friction factor is inversely proportional to the square of the mean velocity. However the 
increase in amplitude has little effect on the mean Nusselt number. Comparing the results of this case with 
the baseline case above, an 18.6% reduction in the axial conduction took place while an increase (by a 
factor of 10.6) in the enthalpy loss occurred. The net effect was about a factor of 6.6 increase in the 
regenerator axial heat loss. The reduction in the axial conduction is attributed to the higher heat flow from 
gas to metal due to the higher instantaneous gas mass flow. This resulted in heating the solid faster in the 
axial direction and thus lower (instantaneous) axial heat conduction. 

Increasing the frequency of the oscillation (by a factor of 3), which also has an effect of increasing the 
velocity amplitude (by a factor of 3), had the effect of decreasing the friction factor upon approaching the 
flow reversal points and increasing the friction factor right after the flow direction has reversed. The mean 
Nusselt number curve was also altered by the increase in frequency of the oscillation. As the flow 
decelerates and approaches the switchover points the mean Nusselt number decreases and when the flow 
accelerates after the direction change, the mean Nusselt number increases. Comparing the results of this 
case with the baseline case, discussed above, we see similar results to the case of increasing the amplitude 
of oscillation. A 14% reduction in the axial conduction took place while an increase (by a factor of 6.9) in 
the enthalpy loss occurred. The net effect was about a factor of 5.8 increase in the regenerator heat loss. 
Again, the reduction in the axial conduction is attributed to the higher heat flow from gas to metal due to 
the higher instantaneous gas flow. 

Increasing the thermal conductivity of the solid material (by changing the solid material from 
stainless steel to nickel) did not have an effect on the friction factor as expected. However, the mean 
Nusselt number showed an overall increase with a more pronounced rise at low Reynolds numbers. It 
should be noted that in this case infinite thermal contact resistance was applied between layers. Therefore, 
this case was compared with the results of the stainless steel material discussed above (not with the base 
case). In this case, a 36.3% increase in the axial conduction took place while a decrease (by 5%) in the 
enthalpy loss occurred. The net effect is about 3.8% increase in the regenerator heat loss. The increase in 
the axial conduction is attributed to the higher thermal conductivity of nickel (compared to stainless 
steel). 

6.10.2 Conclusions: 3-D Straight-Channel-Layer Simulations (Steady and Oscillating Flow) 

For the 3-D straight-channel-layer steady-state simulation, both the friction factor and the mean 
Nusselt numbers depart from agreement with the 2-D simulation values upon entering the second layer. 
That is where the 3-D effects become obvious and they persist as the axial coordinate advances. At the 
entrance of every layer, the forced reorientation of the flow results in small rises of both the friction factor 
and the mean Nusselt number with subsequent decreases as the flow settles into each new layer. Overall 
the plots of the friction factor and the mean Nusselt number tend to flatten out as the flow reaches a fully 
developed condition. 

For the oscillatory-flow simulations of the 3-D straight channel layers, the friction factor shows an 
overall increase compared to the 2-D oscillatory-flow simulation and good agreement with the 
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experimental involute-foil correlation from Gedeon. The mean Nusselt number also shows an overall 
increase compared to the results from the 2-D simulation. It also shows a higher value when compared to 
the correlation from Gedeon at maximum Reynolds number. The shapes of both the friction factor and 
mean Nusselt number curves are similar to the shapes observed in the 2-D simulations. So the effect of 
going from 2-D to 3-D resulted in shifts upwards of both friction factor and mean Nusselt number curves, 
as might be expected since the 2-D simulations do not include the flow perturbations resulting from flow 
around the ends of the foil layer. 

Changing the thermal contact resistance from zero to infinite between the solid layers in the 3-D 
straight-channel layers with oscillating flow had no effect on the friction factor, as expected. However the 
mean Nusselt number increased overall but with a more pronounced increase at the crank angles near 
where the flow switches direction. The same behavior was encountered in the 2-D simulations. Although 
the actual values for the mean Nusselt number are higher in the 3-D simulations, the more expedient 2-D 
simulations capture well the behavior of the mean Nusselt number as the thermal contact resistance is 
changed. 

6.10.3 Conclusions: 3-D Involute-Foil-Layer Simulations (Steady Flow Only) 

Simulations at several Reynolds numbers in the laminar-flow regime were performed and the friction- 
factor and the mean-Nusselt-number results were quantified as functions of their respective dimensionless 
axial coordinates. For Reynolds number 50, comparisons were performed for the 2-D and the 3-D 
straight-channel simulations. The friction-factor variations with Reynolds number were compared to 
experiments and two theoretical correlations. Several conclusions can be drawn: 

As with the 3-D straight-channel layers, the simulations for the 3-D involute-foil layers show 
increases in both friction factor and mean Nusselt number at the geometric transitions between the layers. 
Furthermore, at Reynolds number 50, these increases are similar for the 3-D straight- channel and the 3-D 
involute-foil simulations. From this similarity one can infer that using a simpler grid such as that used for 
the 3-D straight-channel layers can capture to a good extent the steady-state three-dimensional effects of 
the 3-D involute-foil layers for low Reynolds numbers. This is important in terms of the practicality of 
doing CFD simulations. Although more computing power is available now than ever, the researcher still 
must compromise between computing time and accuracy. 

The friction factor matched well the experimental results. That lends credence to the simulations 
performed for the 3-D involute-foil layers. The detailed work that went into constructing the grid, running 
the CFD simulations, and post processing the data provided meaningful results, validated by experiment. 
By their nature, the simulations can in a more expedient way go beyond what the experiments can do and 
provide further predictions for optimization or comparison with other designs. 

The technique of analyzing a repeating unit recursively has also been validated. That allows for 
steady-state simulations of a stack consisting of a large number of layers by using a repeating unit that is 
only two layers thick. That, in turn, not only saves on computation time and resources but also makes the 
simulation of a large stack feasible. 

6.10.4 Future CFD Work 

The CFD involute-foil simulations reported here revealed many insights into how the fluid flow and 
the heat transfer occur in the microfabricated involute-foil-design regenerator. However, the study could 
be usefully extended in two main directions. The first direction consists of varying parameters that result 
in having to modify the grids. Such geometric parameters could be the length of the layer, the hydraulic 
diameter or surface roughness. The second direction consists in changing the boundary conditions and/or 
the material properties. 

As far as boundary conditions are concerned, there are numerous combinations of amplitude and 
frequency that can be attempted based on the conditions encountered by the regenerator in actual Stirling 
engines. Furthermore, inlet profiles that better match the conditions at the entry and the exit of the 
regenerator could be tried. All these new simulations can build upon the present work and save time and 
money, yet provide ever more advanced knowledge. 
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CFD simulations for comparison with the experimental jet penetration study done by UMN have been 
postponed to Phase 111. This study will include a slot jet and a round jet to examine the jet spread angle 
and penetration depth. 


7.0 Structural Analysis of Micro-Fabricated 
Involute-Foil Regenerators (Infinia) 

7.1 Summary 

This section reports the regenerator structural analysis results which helped demonstrate the 
feasibility of the design. The results indicate that the proposed regenerator structure has high axial 
stiffness and the stress level is sensitive to a radial side disturbance. Potential further work is also 
discussed. It should be noted that in these analyses stainless steel was used (the originally planned 
microfabrication material) while the thermal test data (friction factor and Nusselt number) discussed 
earlier was taken using a nickel regenerator. 


7.2 Introduction 

To have high heat storage capacity, regenerators are constructed of high-porosity material that readily 
conducts heat radially and has a high surface area. Most current space-power regenerators are made from 
random fibers, which are somewhat difficult to manufacture in a precisely repeatable manner and are 
susceptible to deformation; these problems can lead to performance losses. The CSU NRA regenerator 
microfabrication contract team proposed a microfabricated involute-foil regenerator to potentially replace 
random-fiber and wire-screen regenerators. Figures 7.1 and 7.2 illustrate the geometry of the annular 
rings and the involute sections of the regenerator. Early analysis showed the potential for significant gains 
in performance efficiency and reductions of manufacturability variability while improving structural 
integrity (Qiu and Augenblick, 2005). 

To ensure that the stiffness and the stress levels meet the design criteria, linear stress analysis was 
carried out on this proposed new regenerator. This section presents the results of finite element analysis of 
the microfabricated, involute-foil regenerator under 44 N (10.0 lb) axial force and 4.4 N (1.0 lb) side 
disturbance force. 



Figure 7.1 . — Partial of the solid model. 
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7.3 Finite Element Analysis 

A finite-element-analysis (FEA) model was created to represent the geometric characteristics and 
three load cases were applied to the finite-element model to examine the stiffness and the stress levels. 

Without using symmetry or periodic-symmetry conditions, all 360° of the geometry was included in 
the finite-element model. The thicknesses of the annular rings and the involute segments were much 
smaller than the other 2 dimensions. To allow the model to be handled by the computer power available, 
the FEA model was simplified as surfaces in 3D space with 4 layers in axial direction. ANSYS shell 
element shell63 was used in the FEA to reduce the size of the model. The regenerator was made from 
stainless steel 3 1 6L and the assumption that the material properties were not sensitive to temperature 
change was used in the analysis. 


7.3.1 Material Properties of Stainless steel 316L Used in FEA 


Young’s modulus 
Poisson’s ratio 
Tensile strength 
Yield strength 


1.9xlO n N/m 2 (2.796xl0 7 psi) 
0.3 

4.97x10 s N/m 2 (7.21 xlO 4 psi) 
1.8x10 s N/m 2 (2.61 xlO 4 psi) 


7.3.2 Geometric Model 


Thicknesses of the involute sections, inside annular rings, and the outside annular ring are 12.7 pm 
(0.0005 in.); 25.4 pm (0.001 in.), and 127 pm (0.005 in.), respectively. Total box volume was about 
270 mm 3 (0.0165 in. 3 ); the mass volume was about 44.5 mm 3 (0.0027163 in. 3 ). The porosity was about 
84%. Individual disk thickness (axial direction) was 250 pm. The FEA was based on the preferred 
stainless-steel material, even though nickel was chosen for convenience in early testing of the 
performance of the involute-foil geometry. 

7.3.3 FEA Model 

ANSYS 3D shell element shell63 with 6 DOF at each node was used. The total number of elements 
was 136422. The total number of nodes was 170220. 

7.3.4 Boundary and Eoading Conditions 
7.3.4.1 Case 1 (Axial Compression) 

It is extremely important that the regenerator be properly fixed within the heater head, as any relative 
movement between the regenerator and the heater head can lead to regenerator structural oscillation and 
failure. Axial compression and a slight press-fit mechanism are currently used by Infinia to stabilize their 
regenerators within the heater head. A 44 N ( 1 0 lb) axial force was uniformly distributed on the top 
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surface to simulate the axial fit and the bottom-face surface was constrained from translation in the axial 
direction. In order to avoid rigid-body motion in the FEA, the minimum constraint condition Ux = Rx = 
Ry = Rz = 0 and Uy = Rx = Ry = Rz = 0 were used on 2 nodes of the inside circle as in the following 
figures. 



Figure 7.3. — FEA model (note 4 layers are modeled). 



Figure 7.4. — Boundary and axial compression loading condition (Case 1). 
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7.3.4.2 Case 2 (Radial Side Load) 

During installation, the regenerator ought to be installed so that it is under axial compression with no 
side loading. To simulate the disturbance from a side load, a 4.45 N (1 lb) side force acting on about 
0.047% of the top layer outside annular ring was added to Case 1 to investigate the side load effect. The 
bottom-face surface was constrained in the axial direction and the inside circle of the bottom face was 
fixed in rotation and translation directions to avoid the rigid-body motion. 

7.3.4.3 Case 3 (Distributed Radial Side Load) 

This load case was created to investigate the stress sensitivity with respect to the side-load acting- 
area. The boundary condition and the load conditions were similar to the Case 2 except that a 4.45 N 
(1 lb) radial side load acted on 10% of the top layer outside annular ring. 

7.3.5 FEA Results for Cases 1, 2 and 3 

The FEA results for Cases 1, 2 and 3 are summarized in table 7.1 and plotted in figures 7.7 through 
7.21. 


TABLE 7. 1— MAXIMUM DISPLACEMENT AND VON MISES STRESS FOR LOAD CASES 1, 2 AND 3 


Load case 

Disp. Ux, 
in. 

{plane of disk} 

Disp. Uy, 
in. 

{plane of disk} 

Disp. Uz, 
in. 

{axial} 

Total disp., 
in. 

von Mises 
stress, 
psi 

Case 1 (axial force) 

0.734e-6 

0.736e-6 

0.646e-6 

0.804e-6 

1732 

Case 2 (rad’l force 1 ) 

0.111e-3 

0.148e-3 

0.289e-5 

0.148e-3 

40624 

Case 3 (rad’l force 2) 

0.462e-4 

0.304e-4 

0.735e-6 

0.462e-4 

6374 



Figure 7.5. — Boundary and radial side-loading condition (Case 2). 
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Figure 7.6. — Boundary and distributed radial side loading condition (Case 3). 
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Figure 7.7. — Deformation Ux (Case 1). 
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Figure 7.8. — Deformation Uy (Case 1). 
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Figure 7.9. — Deformation Uz ( Case 1). 


NASA/CR— 2007-2 15006 


135 


NODAL SOLUTION 


STEP= 
SUB = 
TIHE = 
USUM 
RSYS= 
DHX = 
SMN = 
SHX = 


NODAL 

STEP = 
SUB = 
TIHE = 
SEQV 
DMX = 
SHN = 
SHX = 



AN 

HAR 8 2006 
12:41:14 


Figure 7.10. — Total deformation (Case 1). 
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Figure 7.11. — Von Mises stress (Case 1). 
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Figure 7.12. — Deformation Ux (Case 2). 
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Figure 7.13. — Deformation Uy (Case 2). 
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Figure 7.14. — Deformation Uz (Case 2). 
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Figure 7.15. — Total deformation (Case 2). 
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Figure 7.16. — Von Mises stress (Case 2). 
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Figure 7.17. — Deformation Ux (Case 3). 
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Figure 7.18. — Deformation Uy (Case 3). 
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Figure 7.19. — Deformation Uz (Case 3). 
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Figure 7.20. — Total deformation (Case 3). 
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Figure 7.21. — Von Mises stress (Case 3). 
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7.4 Structural Analysis Summary and Conclusions 

Finite element analysis of the microfabricated involute-foil regenerator shows that the regenerator had 
very high average axial direction stiffness (3.75e7 lb/in.). Without any radial side disturbance, the stress 
level was much lower than the material yielding strength. If the radial side disturbance such as 
misalignment was localized in a small area, as in Case 2, Von Mises stress was beyond the material 
yielding strength and permanent deformation could occur in that area, which may decrease the Stirling 
efficiency. In order to prevent local permanent deformation, the radial side load must be small or the 
disturbance area must be large, as in Case 3. 

In summary, the proposed microfabricated involute-foil regenerator has high axial stiffness. The 
stress level is sensitive to the radial side disturbance, which therefore requires special cautions and 
appropriate processing during installation to prevent lateral permanent deformation. 

7.5 Further Structural Analysis Work 

To enhance the structure of the regenerator in the radial direction and to decrease the likelihood of 
performance degradation should be the main goals of further structural analysis of the microfabricated 
involute-foil regenerator. The stress caused by radial side loading is a function of variables, such as the 
length, the angle, the thickness, the number of the ribs in each involute section, and the number of the 
annular rings. Geometric optimization could be performed to decide the optimum values of the above 
variables for reducing stress. 

8.0 Phase II Conclusions and Recommendations for Future Work 

During Phase II an actual-size regenerator comprised of a stack of 42 disks, 1 9 mm in diameter and 
0.25 mm thick (in the flow direction) — with microscopic involute-shaped flow channels — was 
microfabricated and tested in an oscillating-flow test rig. The geometry resembles an assembly of sections 
of uniformly-spaced parallel plates, except that the plates are curved. The curved sections of plates, or 
“involute foils,” are incoiporated in annular portions of the disks which are separated by concentric rings. 
Two types of disks alternate in the stack, so that the angles between the foils or plates in adjacent disks 
are close to 90°. Each disk was made from electro-plated nickel using the LiGA process. This process 
involved x-radiation of a photoresist through a mask, dissolving portions of the irradiated photoresist, 
then electroplating of nickel on a copper substrate within remaining photoresist channels, etc. This 
regenerator had feature sizes close to those required for an actual Stirling engine, but the overall 
regenerator dimensions were sized for the NASA/Sunpower oscillating-flow regenerator test rig. 
Examination by scanning electron microscope showed the disks were an accurate rendition of the design 
specification, except for a few flaws of types which are expected to be eliminated in the future via 
improvements to the manufacturing process. Testing in the oscillating-flow test rig showed the 
regenerator performed extremely well, producing the highest figures of merit ever recorded for any 
regenerator tested in this rig (since its fabrication about 20 years ago). Other regenerator materials 
recently tested in this rig include random-fiber, wire-screen and etched-foil materials. 

Progress was also made in understanding the detailed fluid dynamics and heat transfer in the 
regenerator by computational fluid dynamic (CFD) analysis at Cleveland State University and large-scale 
testing at the University of Minnesota. In general, the conclusions from the CFD and large-scale testing 
results reinforced those from the actual-size test results and revealed some important details about the 
microscopic flows responsible for the overall regenerator behavior. 

A Phase III effort is now underway to microfabricate a stack of involute-foil disks to form a 
regenerator for testing in a modified Sunpower FTB (Frequency Test Bed) engine. This engine was 
originally designed for a random- fiber regenerator and will not be reoptimized for the new involute-foil 
regenerator — though some modifications will be made to the engine. The Phase III effort includes testing 
of this involute-foil equipped FTB (with hot-end temperature of 650 °C). 
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Beyond this Phase III effort, the microfabrication process needs to be further developed to permit 
microfabrication of higher temperature materials than nickel (Pure nickel is too soft a material for 
practical use even in the 650 °C FTB engine). For example, DOE, NASA and Sunpower are currently 
developing an 850 °C engine for space-power applications. And a potential power/cooling system for 
Venus applications would need regenerator materials capable of -1200 °C. Early Mezzo attempts to 
“EDM” stainless-steel using a LiGA developed EDM tool involved a bum time (dependent on EDM 
machine setting) that was much too large to be practical. Some possible options for further development 
of a micro fabrication process for high-temperature involute-foils are: (1) Optimization of an EDM 
process for high temperature materials that cannot be processed by LiGA only; bum times can be greatly 
reduced by higher-power-EDM-machine settings than originally used in Phase I by Mezzo; but 
“overbum,” i.e., the gaps between the EDM tool and the resulting involute-foil channels, increases with 
higher powers; (2) development of a LiGA only process for some high temperature alloy, or pure metal, 
that would be appropriate for the regenerator application (pure platinum would work but has very high 
conductivity, which would tend to cause larger axial regenerator losses, and it is very expensive), or (3) 
microfabrication of an appropriate ceramic material for high-temperature regenerators (structural 
properties of ceramics, which tend to be brittle, would be a concern). 
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where F L is the configuration factor for the fraction of the total radiation leaving A L that arrives at A(q). 
Similar to the previous case 


f l-\ — j(h -~J H 2 - 4 ) (A.5) 

where 

H = 2 + 4(^ z -0 2 (A.6) 


A.3 Negative Wall Contribution 



q w _ is the radiation flux leaving the wall surface r| < \ that passes through surface A(q). Integrating the 
contributions of differential elements cl A r| it may be written as 

q w- = r ( r i ) 4 F d^ dA ^ 1 ( A - 7 ) 

0 

Substituting for the wall area element 

dA^-4na 2 dr\ (A. 8 ) 

this becomes 

q w- = 4a J 7'(p f Fj^dr \ (A.9) 

0 

/,] - 5 is the configuration factor for the fraction of the total radiation leaving dA n that arrives at A(c). 
From case 30 p. 829 Siegel and Howell 

Vr fr-’V* -fe-i) (A ' 10) 

Vfe - r ,) 2 +1 
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A.4 Positive Wall Contribution 




n 


q w+ is the radiation flux leaving the wall surface r| > \ that passes through surface A(E). Similar to the 
previous case: 


%L 

q w+ = -4 a | r(r)) 4 F+^dr \ (A. 1 1 ) 

% 

F* d n _? is the same configuration factor as the previous case except with c and ?/ switched 

{ r\ + -> -h-a (A.i2) 


A.5 Normalization 

Dividing the 4 heat flux contributions by q mdX = -4 a(T L 4 - T 0 4 ) converts them to normalized radiation 
heat fluxes 

go -{Tq/T l ) 4 F 0 _s 

<7max l-(7o/7i) 4 

U f l-z, 

^max 1 -{To/TlY 

-4 J(A(n)/7/ ) 4 F^dr\ 

Qw- _ 0 

<7 max 1-(To/T l ) 4 

4 J(7 (r,)/7> ) 4 F^dr\ 

*7w+ _ % 

( /niax 1 -{Tq/Tl) 4 


(A. 13) 
(A. 14) 

(A. 15) 

(A. 16) 
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For a linear temperature variation the temperature ratio T(r])/T 0 may be expressed as 


zM 

Tl 


T -^ + 

Tl 


( T \ 

l_Zk JL 

l t l )% l 


(A. 17) 


A.6 Programming 

Custom Delphi program RadiationDownTube.pas performs the above calculations and sums all the 
radiation contributions to produce the results plotted above (Excel does the actual plotting of data). An 
adaptive quadrature routine does the required integrations. The program was tested for two test cases with 
known solutions, (1) the radiation flow in the limit of zero tube length and (2) the radiation flow for a step 
temperature distribution with T = T 0 up to point c, and T = T L beyond. For both cases q/q mm = 1. 
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Appendix B — Regenerator Figure-Of-Merit Degradation With 
Intra-Regenerator Flow-Gap Variations (Gedeon Associates) 


To: Regenerator Research Team 
From: D. Gedeon 
January 7, 2004 

The intra-regenerator flow streaming results for parallel foil regenerators is used to estimate the 
degradation in our regenerator figure of merit as a result of gap variations between sections. 

B.l Background 

We have settled on this formula for calculating the figure of merit for preliminary ranking of 
regenerator matrices 2 


Fm =~ 
f 

For a foil-type regenerator /= 96/Re, Nu = 8.23 and N k = 1. When plotted as a function of Reynolds 
number (assuming Pr = 0.7) the figures of merit for a foil regenerator, as well as a few other regenerators 
of interest, look like the chart below. 

Figures of merit — various matrices 


1 

RePr N k " 

, 4Nu RePr , 



B.2 Effects of Gap Variations 

A foil regenerator is split into two parts. The flow gap in one is increased while decreased in the 
other. This gives rise to DC flow circulations between the two parts and a degradation in overall engine 
efficiency, as previously documented. 3 

Not previously reported was the effect of the gap variations on combined-regenerator pumping loss 
Wf ( Sage output AEfric) and cycle-average enthalpy flow H (Sage output FINeg or FlPos) . Because Wj 
relates to friction factor and H relates to Nusselt number, this information can be used to estimate the 
effective figure of merit as a function of gap variations. Specifically, friction factor is directly 


2 D. Gedeon, Regenerator figures of merit, August 6, 2003, (CSUmicrofabFiguresofMerit.tex) 

3 D. Gedeon, Intra-regenerator flow streaming produced by nonuniform flow channels, Jan 5, 2004 
(CSUmicrofabIntraRegenFlows.doc) 
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proportional to Wf. So the effective friction factor for the gap-perturbed regenerator is greater (or less) 
than the friction factor for the baseline regenerator by the factor 

/ _ WfA + WfB 

7^~{WfA+W fB ) 0 

where the 0 subscript refers to the baseline (equal flow gap) regenerator and subscripts A and B refer to 
the two regenerator parts. Because the Nusselt number is inversely proportional to H , the effective 
Nusselt number for the gap-perturbed regenerator is related to Nusselt number for the baseline regenerator 
by the factor 


Nu _ [h a +H b \ 

Nu 0 H a +H b 

The effective / and Nu thereby computed can be substituted into the formula for figure of merit with 
Reynolds number taken as the baseline regenerator mean value Re = 62. The result is two curves of figure 
of merit degradation as a function of relative gap variation, one for the case where the two regenerator 
parts are in good thermal contact and one for the case of no thermal contact. The latter case is worse 
because of the temperature skewing effect of DC flow which amplifies the DC flow if not suppressed by 
thermal contact between the two parts, as shown in the following plot: 

Parallel Foil Regenerators: 

FM degradation vs flow gap variation 
@ Re = 62 



Delta g / g 0 


The “Delta g / g 0 ” represents the amount by which the gap is higher in part A and simultaneously 
lower in part B, compared to the baseline gap. A gap-variation amplitude, in other words. The curve 
labeled “good thermal connection” corresponds to gap variations between regenerator parts in good 
transverse thermal contact — e.g., between adjacent foil layers or nearby in the same layer. The curve 
labeled “no thermal connection” corresponds to gap variations between distant parts of the regenerator, as 
might occur when all the layers on one side of an annulus are crammed together and spread apart on the 
other side due to some systematic assembly error or misalignment between inner and outer canisters. 

It is seen that in either case it does not take too much gap variation to significantly reduce a foil 
regenerator figure of merit. Increased cyclic enthalpy flow (sum for both parts) is mainly responsible for 
the degradation. Flow frictional loss hardly changes, even decreases a bit. At a gap variation of + 45 % 
the figure of merit is down in the vicinity of random fibers at 90% porosity (F M = 0.14). This is for a 
Stirling engine regenerator. The figure of merit degradation for a cryocooler regenerator would be much 
faster, based on the conclusions of the January 5, 2004, memo. 
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It is reasonable to apply the above plot to the general case of any channel-type regenerator, such as 
the honeycomb type. In that case hydraulic-diameter variation would replace flow-gap variation, but the 
resulting curves should be quite similar. 

So, based on the above, what should we specify for an allowable gap variation (hydraulic diameter 
variation) for a microfabricated regenerator? Maybe a ±10% variation would be reasonable. That should 
keep us above F M = 0.4 for both localized variations and systematic gap variations. If systematic 
variations (over large distances) are not a problem then we might go as high as ± 1 5 %. Maybe a bit more. 
But not too much more lest we wind up in the unenviable position of producing an “improved” 
regenerator that performs as well as random fibers at many times the cost. 
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Appendix C — Potential 6 to 9% Power Increase for a Foil-Type “Microfab” 
Regenerator in the Sunpower ASC Engine (Gedeon Associates) 

To: Regenerator Research Team 
From: D. Gedeon 
March 26, 2004 

Gary Wood has been wondering if the previous estimates for the benefits of a microfabricated 
regenerator (see: August 25, 2003 memo) might have been optimistic considering the regenerator length 
was excessive (133 mm) and foil solid conduction was discounted. Solid conduction multiplier Kmult 
was set to 0.1 as an approximation to the interrupted and convoluted solid conduction path we were 
considering for the lenticular matrix design. 

This memo addresses these concerns. A foil-type is inserted into the canister of the most recent 
Sunpower ASC engine Sage model (D25B3.stl) and the foil spacing is optimized, with length constrained 
first to 70 mm, then 60 mm. The foil is assumed to be 15 pm thick stainless steel with the full solid 
conductivity accounted for in the model. 

After re- optimization the result is a 6.6% increase in PV power for the same heat input for 70 mm 
long foil, compared to the baseline random fiber regenerator. A 5.5% increase for 60 mm long foil. 
Details follow. 


C.l Details 

One of the two Sage files created for this study is named D25B3FoilRegen. It has an optimization 
structure much like that described in a December 23, 2003 memo “Sage Model for ASC Optimization” 
(SunpASCOptimizationFile.doc), except the acceptor length is constrained to 25 mm and the regenerator 
matrix type is foil. Regenerator length is constrained to 70 mm. A second Sage file named 
D25B3FoilRegenL60 is identical except regenerator length is constrained to 60 mm. 

When these Sage files are optimized all of the basic engine dimensions get adjusted to best serve the 
foil regenerator, subject to all the dimensional constraints that have so-far evolved for the ASC engine 
design. The objectives were to see how well a foil regenerator would do compared to the baseline random 
fiber regenerator and get some idea of the trades for reduced regenerator length. Foil gap was optimized 
but foil thickness held constant at 15 pm. 

The table below provides the some key results and dimensions for the two foil regenerator 
optimizations, compared to the baseline random-fiber regenerator. 



Baseline random fiber regenerator 
D25B3 

Foil regenerator 70 mm long 
D25B3FoilRegen 

Foil regenerator 60 mm long 
D25B3FoilRegenL60 

Heat input (W) 

230.0 

230.0 

230.0 

PV power (W) 

111.8 

119.2 

118.0 

Percent increase 

— 

6.6 

5.5 

Foil gap (pm) 

— 

97.2 

92.6 

Canister area (cm 2 ) 

2.83 

2.25 

2.06 


So once again there are significant benefits to foil regenerators, even at reduced length. Note that the 
regenerator canister area went down for the two foil cases resulting in a slightly more compact engine 
with reduced pressure-wall conduction loss. 

For the 60 mm case the foil solid conduction averaged over the regenerator represents about 6.7 W 
out of 230 W heat input. Thi s would be worth another 3% of efficiency (power increase for given heat 
input) if eliminated. For the batch-mode EDM’ed regenerator plates we are considering there would be an 
interrupted solid conduction path allowing us to eliminate some of the 6.7 W conduction loss. So the 
original estimate of a potential benefit on the order of 9% for a microfabricated regenerator continues to 
appear reasonable. 
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Appendix D — EDM Regenerator Disks: Concentric 
Involute Rings (Gedeon Associates) 

To: Microfab Team 
From: D. Gedeon 
April 13,2004 

In which involute foils are packaged in a cylindrical form within concentric rings. 

D.2 Packing a Cylinder 

Thinking ahead to the test regenerator we are planning to build in year 3, we need to come up with a 
matrix that fits within a 19 mm diameter cylindrical form. What should it look like? 

The involute-foil idea is tempting, except that the only way to completely fill a cylinder with an 
involute is with the extreme case of a single foil, wound in spiral fashion. This will not do because it is 
not geometrically similar to the multi-foil involute structure we would use for a thin annular regenerator 
(space power engine) and it results in a long unsupported length of foil, more like a watch spring than a 
regenerator matrix. 

Both problems can be solved by packing the involute elements into concentric rings, each ring 
separated by a thin wall. Each ring may contain a different involute family (different generating circle) so 
they may all be geometrically similar to “thin annular” involutes. The ring walls also serve as support 
points for individual involute elements, thereby increasing their stiffness. 

The illustrations on the following page show two ways one might package involutes into concentric 
rings. In the picture labeled “geometric spacing” the rings are defined by circles with successive 
diameters in the same ratio. The advantage of geometric spacing is that the involutes in all rings have 
about the same angle relative to the cylinder radius and are therefore fluid-dynamically and structurally 
similar insofar as those things matter. The disadvantage is that the ring spacing decreases toward the 
center, resulting in shorter elements (shorter aspect ratios for the flow channels). In the picture labeled 
“uniform spacing” the rings are defined by circles with successive diameters in constant increments. The 
advantage of uniform spacing is that the lengths of the involute elements do not vary as much. The 
disadvantage is that the radius-angle of the involute elements increases toward the center, possibly 
resulting in fluid-dynamic or structural differences between inner and outer rings. 
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Equally spaced involute rings. 


D.2 Dimensions 

The illustrations show cross-section views of regenerator disks created with the Solid-Edge CAD 
software. Only about 14 of the disks are cut into the involute patterns so-as not to bog down the tiny brain 
of Solid-Edge in computing all the little features. Common dimensions for the two regenerator disks are 


Outer diameter 19.05 mm 

Involute channel width (gap) 86 pm 

Wall thickness between involutes «14 pm 

Wall thickness between rings 20 pm 

Disk thickness 500 pm 


The involute channel widths are exact, as these things go, with the wall thickness varying slightly to 
accommodate the approximate circular profiles used instead of true involute profiles. The error so created 
is estimated later on. 


D.3 Central Holes 

In each case there comes a point toward the center of the matrix where continuing the pattern 
becomes absurd and one can either change to another type of foil pattern or just stop and fill the central 
hole with an insulating stuffer. For purposes of the test regenerator it should be acceptable to go with the 
insulating stuffer. For the geometrically-spaced matrix the central hole diameter is 2.6 mm. For the 
equally-spaced matrix it is 5.0 mm. 

D.4 Alternating Sense Involutes 

As illustrated, the angular sense of the involutes alternates in successive rings. This results in a 
herringbone pattern. It would also be possible to maintain the same angular sense in each ring. In either 
case the idea woidd be to flip alternating regenerator disks so the angular sense changes with each layer. 
This way there should be no need for rotational alignment between layers and the solid conduction path 
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between involute wall elements is interrupted. The opportunity for flow distribution between layers is also 
maximized. 

The structural and flow advantages to the herringbone pattern, if any, are not completely clear. On 
possible item of consequence is that any radial flow in a plenum at the discharge end of such a 
herringbone matrix would tend to get mixed because the flow coming from successive involute ring 
would tend to swirl in opposite directions. A good thing, 1 suppose. 

D.5 EDM 

Based on my current understanding of the Mezzo EDM process the illustrated regenerator disks 
should be possible to make. One concession to EDM might be to round the comers where the involute 
elements attach to the circular walls between rings which would result in somewhat greater axial thermal 
conduction loss. 


D.6 Structural Analysis 

In principle one could perform a finite-element stress analysis of a complete involute -ring regenerator 
disk, but 1 have not done that. The important things to understand would be the axial load required to 
buckle a regenerator comprising a stack of such disks and the resistance to deformation of the individual 
involute wall elements. 

The axial buckling load is important if we decide to hold the stack in place by compression. 1 have 
some hope that the axial bucking strength will be adequate because each involute flow channel forms a 
structural cell consisting of side walls integrally connected to end walls. Each structural cell should be 
relatively stiff, compared to foil layers unconnected at the ends. And there are a large number of such 
cells to distribute the load. 

The resistance to deformation of the individual wall elements is important to maintaining uniform 
spacing. If nothing else the walls should be stiff enough so that internal stresses relieved by the EDM 
process do not result in significant spacing variations. 1 am not sure, but it seems that the curved shapes of 
the involute walls will help in this regard. The radius of curvature cannot change too much without 
affecting the separation between endpoints, which is constrained by the end walls (inter-ring walls). Were 
the elements straight (e.g., radial spokes) they could deflect much more for the same end constraints. 
Obviously we need to think more about this. 
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Appendix E — Involute Math (Gedeon Associates) 

Gathered together here are some useful formulas and suggestions for generating involute patterns in 
CAD software. 



Each ring of flow channels is a circular pattern of elemental flow channels. At the center of each 
elemental flow channel is a segment of an involute curve defined by a generating circle and two boundary 
circles, as illustrated in this sketch. 

The generating circle has radius R 0 and the two boundary circles have radii and R 2 . The involute 
segment is the curved arc between R\ and R 2 . By definition it is generated by the mathematical equivalent 
of a string unwinding from the generating circle. Dotted lines indicate the position of the string at the 
beginning and ending of the arc. The exposed lengths of the string at the two endpoints of the involute 
segment are just the arc lengths along the generating circle subtended by angles 9i and 02 or 


n = *001 


and 


r 2 - *002 


r\ and r 2 are also the local radii of curvature of the involute segment at the two endpoints. By drawing 
some right triangles it is easy to conclude that r\ and r 2 are geometrically related to the various circle radii 

by 


n = >?-*« 


and 


A circular pattern of N involute segments can be made by rotating the original involute by angular 
increments 2 n/ N (in radians) about the center of the generating circle. This pattern can also be thought 
of as being made by shortening the string by increments s 0 equal to the circumference of the generating 
circle divided by N, or 
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2tiRq 

N 


The advantage of this point of view is that s 0 is evidently the normal spacing between involute elements. 
Solving the previous equation for N, the number of involute elements spaced by distance So is 


N = 


2kRq 


50 


( 1 ) 


The arc drawn in the above illustration is not actually a true involute curve but rather just an ordinary 
circular arc centered on the generating circle with radius 

r m =(r, +r 2 )/2 (2) 

The normal spacing between the circular arcs is therefore not exactly .v 0 . It varies as 

s® so cos a 

where a is the rotation angle of radius r m relative to its angle at mid segment. This can be seen from the 
following illustration which shows two successive circular “involute” arcs drawn with their centers 
vertically aligned on the generating circle. The reason for the approximately-equals symbol in the 
previous equation is that so is the arc length between two successive center-positions along the generating 
circle, not exactly the cord length. 



E.l Involute Cutouts 

As mentioned, the “involute” arc segments so generated lie at the centers of the flow channels. The 
way the flow channels are generated is to offset the arc segments toward each side by small increments 
(concentric arcs) to create cutout boundaries, as in this illustration of a single flow channel between two 
boundary circles: 
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The “involute” arc segment and boundary circles are colored pink and the cutout boundaries are 
colored blue. In this case the involute arc segment is offset 43 pm toward each side forming a channel of 
width 86 pm. The cutout ends are offset 10 pm from the inter-ring boundary circles, consistent with 
20 pm thick walls between ring sections. The normal distance separating true involute arc segments in the 
circular pattern is 100 pm. With 86 pm channel widths, this leaves 14 pm for the side wall thickness. 
Actually the side walls vary in thickness slightly because of the circular-arc approximation to true 
involute segments and the resulting error in normal spacing that results. That error affects only the side 
wall thickness. The channel width is exact, to within the limits of the CAD software. 

E.2 Spreadsheet Calculations 

Generating the involute channels requires a sequence of boundary circles and a sequence of 
generating circles. In principle the two sequences could be different but it is easiest if the same 
sequence of circles is used for both purposes. In other words, given an arbitrary sequence of circles C„ for 
i = 0. . .M, with diameters D, in increasing order, circle C 0 is the generating circle for the involute arc 
segments between C i and C 2 , circle C 1 is the generating circle for the involute arc segments between C 2 
and C 3 and so forth. Then the only remaining choice is the spacing between the diameters D,. The options 
of geometrical spacing (D,/D, _ t = constant) or arithmetic spacing (D, - D,_ , = constant) have already 
been illustrated and their properties discussed. We might eventually want to look at other options too. 

To automate the process I wrote a spreadsheet that calculates circle sequences and other useful 
information for the 19.05 mm diameter regenerator disks illustrated earlier. Two spreadsheets actually, 
Involute 1 9DiamCyl.xls for geometrically-spaced rings and Involute 1 9DiamCyIEquaISpaced.xls for 
equally-spaced rings. The spreadsheets generate a sequence of circles spanning the desired diameter 
range, then for each concentric ring of the structure calculate the following items, based on the desired 
normal spacing .v (l between involute channels: 

• Number N of elements in the circular pattern (eq. ( 1 ) to nearest integer) 

• Radius of curvature of the approximate-involute arc segment from equation (2) 

• Sweep angle a in radians (©2 — 0] ) 

• Relative spacing error at endpoints 1 - cos a/2 based on equation (3) 

• Spacing error due to roundoff of N 

The last two items are handy for estimating the variation in wall thickness between involute channels. 
For example, the spreadsheet calculated values for the equal-spaced rings are: 
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Circle 

Number 

Radius of 

Sweep angle a, 

Relative spacing error 

Relative spacing error 

diameters D h 

elements in 

involute arcs r m . 

radians 

at endpoints 

due to N roundoff 

mm 

pattern N 

mm 




19.05 






17.05 

472.00 

4.92 

0.24 

0.007 

-0.002 

15.05 

409.00 

4.62 

0.27 

0.009 

-0.002 

13.05 

347.00 

4.29 

0.30 

0.011 

0.000 

11.05 

284.00 

3.94 

0.34 

0.014 

-0.001 

9.05 

221.00 

3.55 

0.40 

0.020 

-0.002 

7.05 

158.00 

3.11 

0.51 

0.033 

-0.004 

5.05 

95.00 

2.60 

0.76 

0.072 

-0.009 


3.05 


Note that the relative spacing error for the inner ring is relatively large (0.072). This times so is the 
approximate amount by which the wall thickness is thinned down near the channel endpoints. For the 
present case of s 0 = 1 00 pm the absolute wall thickness error is about 7 pm. When everything is 
accounted for in Solid Edge the actual wall thickness for the inner ring varies from about 10 pm at the 
ends to 19 pm at mid chord. The wall thickness is much more uniform for subsequent rings. 
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Appendix F — Implications of Regenerator Figure of Merit in 
Actual Stirling Engines (Gedeon Associates) 

Re: NASA/CSU Regenerator Microfabrication Contract NAS3-03124 
From: David Gedeon 
January 27, 2006 


Our regenerator figure of merit measures the heat transfer per unit flow resistance of a regenerator 
matrix. But what does that mean in the context of an actual Stirling engine? (or cooler) The question can 
be answered by imagining a fixed Stirling engine into which regenerators of variable figure of merit F M 
are substituted. It turns out (derived below) that the figure of merit is inversely proportional to the product 
of regenerator pumping loss W p thermal loss Q, and the square of regenerator mean flow area A r : 


F m « 


W p Q t Aj 


For regenerators with the same flow areas the figure of merit is inversely proportional to the product of 
pumping loss and thermal loss. So a high figure of merit will correspond to a low pumping loss, a low 
thermal loss or both. But depending on the relative sizes and importance of the two losses in an actual 
engine the overall benefit to engine efficiency will vary. It is even logically possible that a regenerator 
with higher figure of merit will result in lower actual engine efficiency if it reduces one of the two losses 
that is not very important while allowing the other, that is important, to increase a bit. Or if it reduces the 
engine power density because of a larger void volume. 

F.l Figure of Merit Reformulation 

The figure of merit we have adopted comes from an earlier memo [1]: 


- ■ 


1 


/ 


R e P r 


N k ) 
+ — — 
4Nu Re Pr , 


(1) 


where, 

/ Darcy friction factor 

Nu Nusselt number hd h l k 

N k effective gas conductivity due to thermal dispersion as a fraction of molecular conductivity 
Re Reynolds number put//, /p 

Pr Prandtl number c p \i/ k 

d h hydraulic diameter 

In equation (1) the two terms in the denominator measure the effects of heat transfer and thermal 
dispersion, respectively. Another memo [2] discusses the equivalence between mean-parameter enthalpy 
flows produced by heat transfer and microscopic enthalpy flows produced by thermal dispersion. 

So what does the figure of merit have to do with the losses in a regenerator? To answer that 
question it is convenient to start with the expressions for time average thermal energy transport (enthalpy 
+ dispersion) per unit void flow area q, and pumping power per unit regenerator void volume w r given 
in the 1996 regenerator test rig contractors report[3], appendix C. 
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( 2 ) 


qt = -k 


8T \ Pe 2 Ar 

— i t Nk 

dx \ 4Nu 


and 


Wr = — \fpu 2 \u\) (3) 

where {} stands for time average and Peclet number Pe is shorthand for Re Pr. Multiplying q t by 
regenerator void flow area A r converts it to total thermal energy transport Q t 


Qt = 


1A 8T 

kAf — — 


dx 




(4) 


Multiplying w r by the regenerator void volume A j L converts it to total pumping power W p 


W p = ^{fpui\u\} (5) 

For reasons that are about to become clear it is convenient to express pumping power in terms of FI Re 
by introducing the factor 1 in the form p\u\ dj t /(pRe), which results in 


W p 


A f Lp 2 

2p 



(6) 


One can already see signs of the figure of merit (1) in the preceding equations. It is even clearer by 
writing the figure of merit in the form 


f m - 


r f_ V 

v Rey 


Pr 
Re 2 Pr 2 
4Nu 


+ N k 


(7) 


Ignoring time averages and constants (since we are only interested in proportionalities), representing the 
temperature gradient as AT / L and after a little simplification the figure of merit can be reduced to the 
form 


F m cc 


p 2 c p AT(uAfY 1 

A} WpQt 


( 8 ) 


What does it mean? In the first factor on the right all of the quantities in the numerator are constant for 
any given Stirling machine. 


fixed gas and charge pressure 

fixed hot and cold temperatures 

fixed piston, displacer volumetric flow rates 


c p and p are fixed 

AT is fixed 
uA f is fixed 


NASA/CR— 2007-215006 


164 



So assuming all of the above the figure of merit reduces to 


F m cc 



Or solving for the pumping-work thermal-energy-transport product 


W p Qt*. 


_J 


(9) 


( 10 ) 


Hmm. This interesting result suggests that large regenerator flow areas are good. For two regenerators 
with the same f m the one with the larger flow area will have a lower W p Q t product. 1 suppose that 

makes sense given the dependence of pumping power on velocity to the fourth power in equation (6). 

This may be one reason h\ t does not track overall engine efficiency so closely when comparing 
regenerators of different matrix structures (e.g., parallel plates vs random fibers). 

A cautionary remark in applying equation (10) too literally is that it does not take into account to 
total void volume of the regenerator. The pressure amplitude and therefore power density of the engine 
will vary inversely with the regenerator void volume, though not in direct proportion because there are 
other volumes in the total engine. The implication is that if one measures the W p Q , product relative to the 

square of the engine power then the Aj factor may be somewhat canceled out. But if one is comparing 
two regenerators with the same void volume then this concern vanishes. 
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Appendix G — Random Fibers Correlations with Porosity Dependent 
Parameters — Updated with New 96% Porosity Data (Gedeon Associates) 

Re: NASA/CSU regenerator improvement grant NNC04GA04G 
From: D. Gedeon 
November 1, 2006 

Since the previous memo on this topic 4 additional testing of a 96% porosity random-fiber half-length 
sample “A” produced a lower figure of merit compared to the tests of the full-length 96% porosity 
sample. Recent full-length re-testing suggested that the original full-length data was bad for the 50 bar 
helium test case. 5 

This memo updates the porosity-dependent master correlations for f Nu and N k based on heat transfer 
data for the 96% porosity half-length sample, which is believed to be more reliable. 

1 am planning to go ahead and insert these latest master correlations into the development version of 
Sage because the peak figure of merit for 96% porosity random fibers has dropped significantly from 
about 0.43 down to about 0.28. This will have an effect for far-out regenerator investigations but should 
have minimal impact for analysis of conventional regenerators with porosities on the order of 90%. There 
is a table of relative errors produced by the master correlations near the end of this memo. 

G.l Master Correlations 

The updated correlations are in the same form as before 

/ = — + a-, Re a3 
Re 

Nu = \ + b x ¥e b2 
N k = l + b 3 Pe“ 

Re is the Reynolds number based on hydraulic diameter. Pe = Re Pr is the Peclet number. What has been 
updated is the porosity dependence embedded in the individual coefficients, as follows: (x = |3 /(I - P) , 
where p is porosity) 

=22.7x + 92.3 
a 2 = 0.168x + 4.05 
a 3 =-0.00406x- 0.0759 
b x = (0.00288x + 0.3 10)x 
b 2 =-0.00875jc + 0.631 
*3=1.9 


4 August 19, 2006 memo NASARandomFiberMasterCorrelations.doc) 

5 Details in October 20, 2006 memo NASARandomFiberMysteryLengthDep.doc and November 1, 2006 memo 
N AS ARandomF iberReasonsLengthD ep .doc) . 
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G.2 Figure of Merit 

The overall figure of merit still increases with increasing porosity, but not as much as before. It now 
peaks at a value of about 0.28 at 96% porosity, down from about 0.43 previously. 

0.300 


0.250 


0.200 
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0.100 
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0.000 

10 100 1000 

Reynolds number 

G.3 Correlating Porosity 

The above F u plot and the material in this section is derived from the revised spreadsheet 
RandomFiber Porosity Dependence. xls. 

The data-modeling parameters for the individual regenerator tests are listed in the following table. 

The only changes compared to the table in the August 19, 2006 memo are for the 0.96 porosity Nu and Nk 
parameters, which are now based on tests for half-length sample “A” . 

/parameters Nu, N k parameters 


B 

(3/(l-(3) 

al 

a2 

a3 

bl 

b2 

b3 

0.688 

2.205 

128.8 

3.858 

-0.063 

0.499 

0.635 

3.787 

0.820 

4.556 

248.5 

4.889 

-0.071 

0.945 

0.632 

2.157 

0.850 

5.667 

233.8 

4.15 

-0.082 

1.552 

0.539 

1.113 

0.897 

8.709 

211.2 

5.139 

-0.151 

1.287 

0.600 

1.026 

0.900 

9.000 

321.4 

5.138 

-0.108 

2.323 

0.534 

0.583 

0.930 

13.29 

380.3 

9.906 

-0.195 

7.447 

0.424 

1.983 

0.960 

24.00 

651.5 

6.627 

-0.135 

8.600 

0.461 

2.498 


The individual parameters are plotted below along with their trend-lines. 


Figure of Merit - Master Correlation 
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Friction Factor Parameter 


0.000 5.000 10.000 15.000 20.000 25.000 30.000 

porosity/(1-porosity) 





Nu, Nk Parameter 



porosity/(1 -porosity) 


Friction Factor Parameter 


Nu, Nk Parameter 



porosity/(1 -porosity) 


4 
3 

a 2 

i 

0 

0.000 5.000 10.000 15.000 20.000 25.000 30.000 

porosity/(1 -porosity) 


♦ y = -5.24E-03x + 1 .93E+00 

or v = 1.9 about as good 

♦ 


♦ ♦ 
♦ 


The previous quadratic curve fit (trend-line) for the b\ heat-transfer parameter (Nusselt number 
coefficient) is no longer obviously necessary but it does reduce relative errors for the individual 
correlations, compared to a linear trend line, and also results in a more uniform progression of increasing 
figure of merit with porosity. The y-intercept remains zero so that the implied heat-transfer coefficient 
does not diverge at zero porosity (h oc Nu (1 - p)/p) as argued in the August 19, 2006 memo. 

The following table gives the updated table of correlated vs individual ratios for various quantities of 
interest, averaged over the Reynolds number range from 10 to 1000. Values for individual correlations are 
indicated by zero subscripts. F M is the overall figure of merit. 
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Master correlation relative errors 
Averaged over Re = 10 to 1000 

p 

f/fo 

Nu/Nuo 

Nk/Nu, 

Fm/F mo 

0.688 

1.05 

1.22 

0.47 

1.35 

0.820 

0.83 

1.28 

0.75 

1.52 

0.850 

1.03 

1.40 

1.95 

1.09 

0.897 

1.35 

1.80 

1.47 

0.99 

0.900 

0.99 

1.39 

3.09 

0.99 

0.930 

0.97 

0.93 

1.40 

0.94 

0.960 

0.98 

0.89 

0.66 

1.12 


Compared to the table in the August 19, 2006 memo (including 0.93 porosity data) the figure-of-merit 
ratios in the range of porosities 0.688 to 0.85 are substantially higher. There is not much change at 
porosity 0.90 and a bit more pessimism at 0.93. The correlation at porosity 0.96 is based on new data. 

G.4 Data Files 

The derived data files used for parameter modeling were: 


Sample 

Tested 

Porosity 

DP file 

HX file 

2 mil Brunswick 

1992-93 

0.688 

P06-29Scaled 

H 1 1 -2 1 Scaled 

1 mil Brunswick 

1992-93 

0.82 

PI l-04Scaled 

HI l-18Scaled 

30 pm Bekaert 

2006 

0.85 

DP 85Porosity 

HX 85PorosityTrunc 

12 pm Bekaert 

2003 

0.897 

BekD 1 2P90DPRetestScaled 

BekD 12P90HXScaled 

30 pm Bekaert 

2006 

0.90 

DP 90Porosity 

HX 90PorosityTrunc 

30 pm Bekaert 

2006 

0.93 

DP 93Porosity 

HX 93PorosityTrunc 

30 pm Bekaert 

2006 

0.96 

DP 96Porosity 

HX 96PorosityHalfA 
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